Primer Corte

Unidad 1: Que es el BA

El Business Analytics es el uso de los datos para generar insights que pueden ayudar a la toma de decisiones en las empresas. Para poder realizar este proceso hay que partir de una pregunta de negocio, como “Como puedo incrementar las ventas del cereal?” o “Que producto le interesa más a mis clientes?”, estas preguntas guiaran el proceso de Business Analytics por cuatro “etapas”, que se denominan: tipos de analítica, Estos son:

Descriptiva: que esta pasando (business intelligence)
Diagnóstica: por que esta pasando
Predictiva: que va a pasar
Prescriptiva: que debo hacer

Dentro de estos tipos de analítica se realizan unas designadas tareas las cuales nos darán los resultados para crear los insights o las conclusiones, estas tareas de la analítica son:

Clasificación: predecir la categoría de un individuo
Regresión: predecir un valor continuo
Detección de Excepciones: encontrar anomalías
Clustering: agrupar individuos por elementos similares
Asociaciones: encontrar reglas de coocurrencia
Pronóstico: regresión con un componente temporal
Resumen:simplificar la representación de información
Visualización: facilitar la comprensión de los datos por medio de visualizaciones

Y estas tareas entran dentro de los tipos de la siguiente manera:

Ahora bien, todo esto no lo hace solo una persona, de hecho, hay un grupo entero de personas que se encargan de una parte del proceso del business analytics (y business intelligence). Pero antes de mencionarlos, cual es este proceso? Bueno el proceso del BA es el siguiente:

Recolección de datos y almacenado
Preparación de los datos
Exploración y visualización
Modelado, experimentación, predicción, etc.
Comunicación de resultados
Decisiones

Con los pasos del BA ya definidos, podemos ahora si definir los roles de la analítica y cual es el trabajo que desempeñan. Antes de empezar a mencionar al equipo de analítica, podemos empezar con el analista de datos, que es el encargado de los trabajo de BI, es decir, se queda únicamente en las tareas descriptivas y usa principalmente herramientas como PowerBI, Excel y SQL; se puede ser analista de datos con un curso. Ahora los tres personajes principales en los roles de la analítica son:

Ingeniero de Datos: encargado de preparar los datos para que el científico de datos los analice. Trabaja en el primer paso de la analítica y utiliza SQL, Java, Scala y/o Python; estos ingenieros generalmente son ingenieros de sistemas.
Científico de Datos: encargado de estimar modelos estadísticos y/o de IA, es decir, el encargado de realizar las tareas mencionadas anteriormente. Para esto, utiliza principalmente SQL, Python y/o R y se enfoca en el segundo al cuarto paso de la analítica, incluso puede realizar el quinto pero si no es tan autista como lo son la mayoría; estos científicos son economistas con una maestría.
Analytics Translator: Como ya dije, no siempre el científico de datos va a poder adaptar lo que aplico al lenguaje de alguien que no sabe lo que hizo, por lo que existe este rol, que tiene un conocimiento de la ciencia de datos y de los negocios, lo que lo ayuda a traducir los modelos a insights.

Unidad 2: Introducción a R y Tidyverse

Ahora si, dejando de un lado un poco la teoría, para poder interpretar el BA, y más importante, para poderlo realizar, es necesario realizar código para poder hacer estos modelos y hallazgos. Es por esto que empezaremos a analizar que es R y que podemos hacer con los primeros paquetes, dplyr y ggplot2, que provienen de Tidyverse, para poder hacer análisis básicos de datos.

Para este análisis comenzaremos cargando nuestras librerías y la base que usaremos, que en este caso sera gapminder:

require(pacman)
p_load(tidyverse, gapminder)

data("gapminder") #Si les sale <promise> le dan click y les carga normal

Con esto ya tenemos el paquete que deseamos utilizar, pero antes de poder empezar a aprender como usar Tidyverse, es importante aprender algunas funciones básicas de R, como lo son la creación de objetos, de listas, de matrices, de data frames y secuencias:

objeto <- 10
lista <- c(1, 2, 3)
matriz <- matrix(c(1, 2, 3, 4), nrow = 2, ncol = 2)
data.frame <- data.frame(data1 = c(1, 2, 3, 4), 
                         data2 = c(5, 6, 7, 8))
secuencia1 <- c(1:5)
secuencia2 <- seq(1, 5, by = 2) #Hace lo mismo que x:y solo que se puede cambiar conteo de la secuencia

Ahora, ya que sabemos como crear listas y objetos, podemos aprender que hacer con estos, como ya dijimos que BA se basa en modelos estadísticos, se pueden usar funciones que hagan, medias, medianas y otras estadísticas descriptivas:

prueba <- rnorm(15) #Creamos una lista de 15 numeros distribuidos normalmente

mean(prueba)

## [1] -0.03104879

median(prueba)

## [1] 0.05600307

sd(prueba)

## [1] 1.297628

summary(prueba)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.60196 -0.75319  0.05600 -0.03105  0.78052  1.88061

Estas unas de las funciones que podemos utilizar para hacer estadística descriptiva básica en R, ya para hacer operaciones como suma, resta, división, etc; se hace exactamente igual a como se hace en Excel. Ahora con eso ya fuera del camino podemos empezar a aprender como utilizar Tidyverse, especialmente dplyr y ggplot.

Antes de aprender como usar estos dos, es necesario saber que tanto en ggplot como en dplyr se pueden hacer cadenas de código de una vez en lugar de correr las cosas una por una, es por esto que si se esta realizando alguna tarea con dplyr, utilizamos el conector %>% y si es con ggplot usamos +. A continuación mostrare las funciones más importantes de dplyr.

gapminder %>%
  filter(continent == "Americas") #Nos da solo las observaciones que son del continente de las Americas

## # A tibble: 300 × 6
##    country   continent  year lifeExp      pop gdpPercap
##    <fct>     <fct>     <int>   <dbl>    <int>     <dbl>
##  1 Argentina Americas   1952    62.5 17876956     5911.
##  2 Argentina Americas   1957    64.4 19610538     6857.
##  3 Argentina Americas   1962    65.1 21283783     7133.
##  4 Argentina Americas   1967    65.6 22934225     8053.
##  5 Argentina Americas   1972    67.1 24779799     9443.
##  6 Argentina Americas   1977    68.5 26983828    10079.
##  7 Argentina Americas   1982    69.9 29341374     8998.
##  8 Argentina Americas   1987    70.8 31620918     9140.
##  9 Argentina Americas   1992    71.9 33958947     9308.
## 10 Argentina Americas   1997    73.3 36203463    10967.
## # ℹ 290 more rows

gapminder %>%
  select(continent) #Seleccionamos la variable continente

## # A tibble: 1,704 × 1
##    continent
##    <fct>    
##  1 Asia     
##  2 Asia     
##  3 Asia     
##  4 Asia     
##  5 Asia     
##  6 Asia     
##  7 Asia     
##  8 Asia     
##  9 Asia     
## 10 Asia     
## # ℹ 1,694 more rows

Es importante saber que para estas dos funciones se pueden hacer las dos a la vez y varias condiciones, por ejemplo:

gapminder %>%
  filter(continent == "Americas", year == 2007) %>%
  select(country, gdpPercap)

## # A tibble: 25 × 2
##    country            gdpPercap
##    <fct>                  <dbl>
##  1 Argentina             12779.
##  2 Bolivia                3822.
##  3 Brazil                 9066.
##  4 Canada                36319.
##  5 Chile                 13172.
##  6 Colombia               7007.
##  7 Costa Rica             9645.
##  8 Cuba                   8948.
##  9 Dominican Republic     6025.
## 10 Ecuador                6873.
## # ℹ 15 more rows

Esas son las funciones que se utilizan para el manejo de bases de datos, pero también existen otras dos que se utilizan para procedimientos, una se utiliza para realizar procedimientos sobre las observaciones (summarize) y la otra se utiliza para crear variables:

gapminder %>%
  summarize(prom = mean(gdpPercap)) #Solo nos da el resultado

## # A tibble: 1 × 1
##    prom
##   <dbl>
## 1 7215.

gapminder %>%
  mutate(GDP = gdpPercap*pop/1000000) %>%
  mutate(xd = "Si") %>%
  select(GDP, xd) #Es necesario selecionar la variable para tener solo el resultado

## # A tibble: 1,704 × 2
##       GDP xd   
##     <dbl> <chr>
##  1  6567. Si   
##  2  7585. Si   
##  3  8759. Si   
##  4  9648. Si   
##  5  9679. Si   
##  6 11698. Si   
##  7 12599. Si   
##  8 11821. Si   
##  9 10596. Si   
## 10 14122. Si   
## # ℹ 1,694 more rows

Ahora, algo para complementar la función de summarize es el uso de otra función llamada group_by que como su nombre lo indica es para agrupar observaciones iguales, esto nos es útil porque si notaron el uso de summarize previo, nos da el promedio de todos los gdpPercap de todos los continentes en todos los años, por lo que se puede seleccionar un año y tener el promedio por continente de la siguiente manera:

gapminder %>%
  filter(year == 2007) %>%
  group_by(continent) %>%
  summarize(prom = mean(gdpPercap))

## # A tibble: 5 × 2
##   continent   prom
##   <fct>      <dbl>
## 1 Africa     3089.
## 2 Americas  11003.
## 3 Asia      12473.
## 4 Europe    25054.
## 5 Oceania   29810.

Ahora, ya que cubrimos las bases de dplyr, podemos aprender a usar ggplot, y es que si nos damos cuenta, tener todas esas tablas y números ahí no nos ayudan mucho como a la interpretación, es por eso que con el uso de ggplot2 se pueden realizar visualizaciones para hacer el análisis de los datos más fácil.

Hay 9 gráficos fundamentales para la tarea de visualizaciones en ggplot, estos son:

Histograma: para 1 variable cuantitativa (numérica/integer)
Boxplot: para 1 variable cuantitativa pero se puede usar una variable cualitativa (character/factor) para hacer varios en un plot
Barras: para 1 variable cualtitativa
Columnas: gráfico de barras pero con una variable cuantitativa
Barras Aplicadas/Agrupadas: gráfico de barras con varias observaciones para una variable cuantitativa que se puede hacer en grupos
Puntos/Jitter: para dos variables cuantitativas, pero se puede usar una cualitativa para el color o para el tamaño, o para las dos
Linea: para dos variables cuantitativas, donde una puede ser tiempo
Smooth: para hacer regresiones estadísticas (se hace junto a un jitter)

Para realizar una visualización se deben tener en cuenta las siguientes capas (de abajo hacia arriba):

Donde:

Datos: los datos que se utilizaran para el gráfico
Aesthetics: las variables que se utilizaran para el gráfico
Geometría: el tipo de gráfico
Facets: generar subgráficos en el plot
Estadísticas: añadir información estadística a un gráfico (como geom_smooth)
Escalas: definir la escapa de los ejes
Coordenadas: definir como se van a definir los ejes y el espacio donde estarán
Tema: la apariencia final del gráfico

A continuación se hará un ejemplo de cada uno, por lo que crearemos 3 data frames para hacer los gráficos:

d1 <- filter(gapminder, year == 2007)

d2 <- gapminder %>% 
  group_by(year, continent)

d3 <- gapminder %>% 
      group_by(year, continent) %>% 
      summarise( Poblacion = sum(pop)) %>% 
      group_by(year) %>% 
      mutate(Prop_Poblacion = 100 * Poblacion / sum(Poblacion))

Ahora si, comenzamos con cada uno, en el orden ya establecido:

ggplot(d1, aes(x = pop)) +
  geom_histogram(col = "black", fill = "purple") +
  labs(title = "Histograma de la Población de los Paises en 2007",
       x = "Población", y = "Total") +
  theme_bw()

ggplot(d1, aes(x = continent,y = pop)) +
  geom_boxplot(col = "black", fill = "purple", alpha = 0.7) +
  labs(title = "Boxplot de la Población por Continentes en 2007",
       x = "Continente", y = "Población") +
  theme_bw()

ggplot(d1, aes(x = continent)) +
  geom_bar(col = "black", fill = "purple") +
  geom_label(stat = "count", aes(label = ..count..)) + #Un extra, que es ponerle el numero a las barras
  labs(title = "Total de Paises por Continente", subtitle = "En el año 2007",
       x = "Continente", y = "Total", caption = "Fuente: Elaboración Propia") +
  theme_bw()

d1 %>%
  filter(continent == "Oceania") %>%
  ggplot(aes(x = country, y = pop)) +
  geom_col(col = "black", fill = "purple") +
  geom_label(aes(label = round(pop/10000, 2)), hjust = 0.75) +
  theme_bw() +
  coord_flip()

ggplot(d2, aes(x = as.character(year), y = pop, fill = continent))+
  geom_bar(position = "dodge", stat = "identity") +
  labs(y = "Total de Población", 
       x = "Año", fill = "Continente") +
  theme_bw()

ggplot(d3, aes(x = as.character(year), 
               y = Prop_Poblacion, fill = continent)) +
  geom_bar(position = "stack", stat = "identity") +
  labs(y = "% de la población mundial", 
       x = "Año", fill = "Continente") +
  theme_bw()

ggplot(d1, aes(x = gdpPercap,y = pop)) +
  geom_point(col = "purple") +
  theme_bw()

ggplot(d1, aes(x = gdpPercap,y = pop)) +
  geom_jitter(col = "purple", alpha = 0.7) +
  theme_bw()

gapminder %>%
  filter(country == "Colombia") %>%
  ggplot(aes(x = year, y = pop)) +
  geom_line(col = "purple", lwd = 1) +
  theme_bw()

ggplot(d1, aes(x = gdpPercap,y = pop)) +
  geom_jitter(col = "purple", alpha = 0.7) +
  geom_smooth(data = d1, aes(x = gdpPercap, y = pop), method = "lm", 
              col = "black") +
  theme_bw()

Con todo eso, tenemos todos los gráficos básicos que se pueden realizar, de igual manera es posible mejorarlos utilizando facets que hacen que hallan varios gráficos en un plot, lo que puede ser util a la hora de explicar algunos datos. Por ejemplo:

ggplot(d1, aes(x = gdpPercap,y = pop)) +
  geom_jitter(col = "purple", alpha = 0.7) +
  facet_wrap(continent ~., scales = "free") +
  theme_bw()

Otra manera de hacer es usando un grid.arrange del paquete gridExtra, de la siguiente manera:

#install.packages("gridExtra")
library(gridExtra)

p1 <- ggplot(d1, aes(x = gdpPercap,y = pop)) +
  geom_jitter(col = "purple", alpha = 0.7) +
  theme_bw()


p2 <- ggplot(d1, aes(x = gdpPercap,y = pop)) +
  geom_smooth(method = "lm", col = "black") +
  theme_bw()

grid.arrange(p1, p2)

Es necesario saber que no siempre un gráfico va a ser perfecto, e incluso, muchas veces es posible que el proceso de una realizar una buena visualización sea más largo de que uno espera, es por esto que es importante saber unos tips para mejorar las visualizaciones. Estos tips son:

Organizar los datos: muchas veces es posible que un gráfico se vea mal porque hay muchas observaciones entonces no se puede leer los ejes. En estos casos se pueden voltear los ejes o incluso organizarlos en orden alfabético
Evitar colores que no comunican: Si hay muchos colores y estos no dicen nada, lo mejor es simplemente elegir un color para todo el gráfico o incluso resaltar una observación utilizando, en el caso de gapminder:
```
gghighlight(country == 'Colombia', use_direct_label = F)
```
Evitar gráficos spaghetti: son gráficos que tienen lineas demasiado juntas, donde literal no se puede interpretar nada, en estos casos se puede hacer lo mismo que se mostró con el gghighlight. Estos gráficos se ven así:
No utilizar gráficos torta: no comunican nada y se puede explicar mejor con un gráfico de barras
Gráficos interactivos: con la función plot_ly del paquete plotly se pueden hacer gráficos interactivos donde podes ver que es la observación, podes desactivarla o activarla, pero cuidado con exceder su uso.

Al finalizar cada unidad vamos a limpiar el workspace para que luego no quede demasiado lleno:

rm(list = ls())

Unidad 3: Otras Visualizaciones

Aparte de los gráficos de ggplot que ya vimos, existen otros gráficos que se pueden realizar con otros paquetes, estos gráficos son:

Treemap

Estos gráficos sirven para resumir 3 o 4 variables , donde dos son variables cuantitativas y una o dos son cualitativas. Estos gráficos se ven así:

#Libreria
p_load(tidyverse, treemap, treemapify)

ggplot(gapminder.2007, aes(fill = lifeExp, area = pop, label = country)) + geom_treemap() +
geom_treemap_text(colour = "white", place = "centre")

Esto es en el caso de usar 2 variables cuantitativas y 1 cualitativa, pero si deseamos usar dos variables de cada tipo, el gráfico se vería de esta manera:

ggplot(gapminder.2007, aes(fill = lifeExp, area = pop, label = country, 
                           subgroup = continent)) + geom_treemap() +
  geom_treemap_subgroup_border() + 
  geom_treemap_subgroup_text(place = "centre", grow = T, alpha = 0.5, 
                             colour ="black", fontface = "italic", min.size = 0) +
  geom_treemap_text(colour = "white", place = "centre")

Gráfico de Embudo

El gráfico de embudo se utiliza para mostrar los datos en diferentes etapas de un proceso de negocio, estos gráficos se pueden utilizar para mostrar procesos que tiene falencias o los distintos pasos que da un usuario al entrar a un sitio web. A continuación se muestran dos aplicaciones de un embudo:

#Libreria
p_load(plotly)

fig <- plot_ly() 
fig <- fig %>%
  add_trace(type = "funnel",
            y = datos.funnel$Fase ,
            x = datos.funnel$Clientes,
             textposition = "inside",
            textinfo = "value+percent initial",
            opacity = 0.65)
fig <- fig %>%
  layout(yaxis = list(categoryarray = datos.funnel$Fase))

fig

En este caso solo usamos una variable, pero es posible incluir más variables en el análisis, lo que hace que la visualización se vea de esta manera:

ibrary(plotly)
fig2 <- plot_ly(
    type = "funnel",
    name = 'Prensa',
    y = datos.funnel2$Etapa,
    x = datos.funnel2$Prensa,
    textinfo = "value+percent initial")

fig2 <- fig2 %>%
  layout(yaxis = list(categoryarray = datos.funnel2$Etapa))

fig2 <- fig2 %>%
  add_trace(type = "funnel",
    name = 'Radio',
    y = datos.funnel2$Etapa,
    x = datos.funnel2$Radio,
    textinfo = "value+percent initial",
    textposition = "inside"
    )

fig2 <- fig2 %>%
  add_trace(type = "funnel",
    name = 'Buscadores',
    y = datos.funnel2$Etapa,
    x = datos.funnel2$Buscadores,
    textinfo = "value+percent initial",
    textposition = "inside"
    )

fig2 <- fig2 %>%
  add_trace(type = "funnel",
    name = 'Links',
    y = datos.funnel2$Etapa,
    x = datos.funnel2$Links,
    textinfo = "value+percent initial",
    textposition = "inside"
    )

fig2 <- fig2 %>%
  add_trace(type = "funnel",
    name = 'Directa',
    y = datos.funnel2$Etapa,
    x = datos.funnel2$Directa,
    textinfo = "value+percent initial",
    textposition = "inside"
    )
fig2

Diagrama de Sankey

El diagrama de Sankey se utiliza para ver el cambio de estados de unos individuos entre periodos, por ejemplo, el gráfico que voy a mostrar ahora es un gráfico que representa la migración de las personas de una región (de la izquierda) a otra (a la derecha):

p_load(tidyverse, viridis, patchwork, hrbrthemes, circlize, networkD3)

nodes <- data.frame(name=c(as.character(data_long$source), as.character(data_long$target)) %>% unique())
 
data_long$IDsource=match(data_long$source, nodes$name)-1 
data_long$IDtarget=match(data_long$target, nodes$name)-1

ColourScal ='d3.scaleOrdinal() .range(["#FDE725FF","#B4DE2CFF","#6DCD59FF","#35B779FF","#1F9E89FF","#26828EFF","#31688EFF","#3E4A89FF","#482878FF","#440154FF"])'

sankeyNetwork(Links = data_long, Nodes = nodes,
              Source = "IDsource", Target = "IDtarget",
              Value = "value", NodeID = "name", 
              sinksRight=FALSE, colourScale=ColourScal, nodeWidth=40, 
              fontSize=13, nodePadding=20)

Gráfico de Calendario (Calendar Plot)

Estos gráficos se utilizan en las series de tiempo para analizar regularidades (patrones) en los días de la semana, y sacar insights que no se conseguirían en gráficos normales, como hallar un patrón de comportamiento todos los lunes. Estos gráficos se verían de esta manera:

p_load(tidyquant, openair, RColorBrewer)

retornos.diario <- Return.calculate(PORTAFOLIO2, method = "log")*100 

date = index(retornos.diario)
RETORNOS.DIARIO.DF = data.frame(date, retornos.diario)

colores <- brewer.pal(9, "PuBu")
calendarPlot(RETORNOS.DIARIO.DF, pollutant = "GRUPOSURA", year = 2017, cols= colores)

Segundo Corte

Unidad 4: Limpieza y Exploración de Datos

Como se explico al comienzo, una de las partes en el proceso del BA es la limpieza y la exploración de datos, de hecho, este paso puede tomar hasta el 80% del proceso entero. Este proceso consiste en ver la estadística descriptiva de cada una de las variables y de buscar errores en la base como lo pueden ser duplicados, datos fuera de rango y outliers. Para esta unidad necesitaremos los siguientes paquetes:

require(pacman)
p_load(tidyverse, forcats, reshape2, funModeling, Hmisc, finalfit, rio)

diamonds <- import("bases/diamonds.csv") #Import es una función de Rio, pueden usar read.csv

Algo fundamental a la hora de la exploración de datos, es que siempre es necesario ver los datos en caso de que halla una variable que no aporte nada, poderla eliminar, de igual manera, es importante ver las clases de las variables para ver que si estén bien categorizadas:

head(diamonds, 5)

##   V1 carat     cut color clarity depth table price    x    y    z
## 1  1  0.23   Ideal     E     SI2  61.5    55   326 3.95 3.98 2.43
## 2  2  0.21 Premium     E     SI1  59.8    61   326 3.89 3.84 2.31
## 3  3  0.23    good     E     VS1  56.9    65   327 4.05 4.07 2.31
## 4  4  0.29 Premium     I     VS2  62.4    58   334  4.2 4.23 2.63
## 5  5  0.31    Good     J     SI2  63.3    58   335 4.34 4.35 2.75

sapply(diamonds, class)

##          V1       carat         cut       color     clarity       depth 
##   "integer"   "numeric" "character" "character" "character"   "numeric" 
##       table       price           x           y           z 
##   "numeric"   "integer" "character"   "numeric"   "numeric"

Viendo esto, nos damos cuenta de que hay una variable que no es necesaria (V1) porque solo nos dice el numero del la observación, pero para eso ya tenemos la filas numeradas, de igual manera, tenemos que cut, color, clarity y x estan mal categorizadas, por lo que vamos a hacer eso primero:

diamonds <- diamonds[-1] #se puede hacer con select pero esto es más rapido

diamonds <- diamonds %>%
    mutate(x = as.numeric(x)) %>%
  mutate_if(is.character, as.factor)

Con esto hecho, ya podemos empezar la exploración preliminar de los datos y para esto debemos primero separar las variables entre cuantitativas y cualitativas. Afortunadamente con el paquete de funModeling puede hacer esto automáticamente:

var <- data_integrity(diamonds)

v_cuanti <- var$results$vars_num
v_cuanti

## [1] "carat" "depth" "table" "price" "x"     "y"     "z"

v_cuali <- var$results$vars_cat
v_cuali

## [1] "cut"     "color"   "clarity"

Comenzamos primero analizando las variables cuantitativas, donde para este tipo el proceso a realizar es:

Análisis Estadístico
Análisis Gráfico
Valores Atípicos

Con eso dicho, para comenzar el análisis estadístico, vemos primero si hay datos que son demasiado fuera de lo normal, es decir, más atípicos que los atípicos, para esto primero usamos la función describe del paquete Hmisc que nos mostrara esto:

describe(diamonds[v_cuanti]) #El [] le dice a R que debe ser algo especifico

## diamonds[v_cuanti] 
## 
##  7  Variables      53940  Observations
## --------------------------------------------------------------------------------
## carat 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##    53940        0      274    0.999   0.7985     0.75   0.5132     0.30 
##      .10      .25      .50      .75      .90      .95 
##     0.31     0.40     0.70     1.04     1.51     1.70 
## 
## lowest : 0.2  0.21 0.22 0.23 0.24, highest: 4.01 4.13 4.5  5.01 29  
## --------------------------------------------------------------------------------
## depth 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##    53940        0      185    0.999    61.75     61.8    1.517     59.3 
##      .10      .25      .50      .75      .90      .95 
##     60.0     61.0     61.8     62.5     63.3     63.8 
## 
## lowest : 1    43   44   50.8 51  , highest: 72.2 72.9 73.6 78.2 79  
## --------------------------------------------------------------------------------
## table 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##    53940        0      127     0.98    57.46     57.5    2.448       54 
##      .10      .25      .50      .75      .90      .95 
##       55       56       57       59       60       61 
## 
## lowest : 43   44   49   50   50.1, highest: 71   73   76   79   95  
## --------------------------------------------------------------------------------
## price 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##    53940        0    11603        1     3933     3180     4012      544 
##      .10      .25      .50      .75      .90      .95 
##      646      950     2401     5324     9821    13107 
## 
## lowest :    23   326   327   334   335, highest: 18803 18804 18806 18818 18823
## --------------------------------------------------------------------------------
## x 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##    53939        1      554        1    5.731    5.685    1.276     4.29 
##      .10      .25      .50      .75      .90      .95 
##     4.36     4.71     5.70     6.54     7.31     7.66 
## 
## lowest : 0     3.73  3.74  3.76  3.77 , highest: 10.01 10.02 10.14 10.23 10.74
## --------------------------------------------------------------------------------
## y 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##    53940        0      554        1    5.736     5.69    1.273     4.30 
##      .10      .25      .50      .75      .90      .95 
##     4.36     4.72     5.71     6.54     7.30     7.65 
## 
## lowest : 0     2     3.68  3.71  3.72 , highest: 10.16 10.54 31.8  58.9  90   
## --------------------------------------------------------------------------------
## z 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##    53940        0      376        1    3.539     3.51   0.7903     2.65 
##      .10      .25      .50      .75      .90      .95 
##     2.69     2.91     3.53     4.04     4.52     4.73 
## 
## lowest : 0    1.07 1.41 1.53 2.06, highest: 6.72 6.98 8.06 10   31.8
## --------------------------------------------------------------------------------

Ahí nos dice que hay dos atípicos, y son un 29 en carat y un 90 en y, esto se puede concluir, por el que valor anterior a el es un 5.01 y un 58.9 respectivamente. Antes de eliminarla veamos la información de los atípicos:

diamonds[which.max(diamonds$carat),]

##      carat       cut color clarity depth table price    x    y    z
## 1379    29 Very Good     G     SI2  62.6    56  2967 6.01 6.07 3.78

diamonds[which.max(diamonds$y),]

##     carat  cut color clarity depth table price    x  y    z
## 802  1.01 Fair     E     SI2  67.6    57  2862 6.21 90 4.18

En este caso, vamos a asumir que es un error de digitación, por lo que diremos que ese 29 es un 2.9 y ese 90 es un 9.0:

diamonds$carat[which.max(diamonds$carat)] <- 2.9
diamonds$y[which.max(diamonds$y)] <- 9.0

Ahora utilizaremos unas funciones de funModeling que nos darán la estadística descriptiva de las variables, estas son profiling_num que nos dará esta estadística descriptiva por variable, pero si se desea, se puede hacer de una sola variable. La otra función va a ser df_status que nos dará el total de 0s, porcentaje de 0s, total de NAs, porcentaje de NAs, total de valores infinitos, porcentaje de valores infinitos, tipo de variable y total de valores únicos en la variable:

profiling_num(diamonds)

##   variable         mean      std_dev variation_coef   p_01   p_05   p_25
## 1    carat    0.7979776    0.4740976     0.59412393   0.24   0.30   0.40
## 2    depth   61.7482833    1.4563039     0.02358452  57.90  59.30  61.00
## 3    table   57.4571839    2.2344906     0.03888966  53.00  54.00  56.00
## 4    price 3932.7898776 3989.4487309     1.01440678 429.00 544.00 950.00
## 5        x    5.7311841    1.1217537     0.19572809   4.02   4.29   4.71
## 6        y    5.7344577    1.1425660     0.19924569   4.04   4.30   4.72
## 7        z    3.5388730    0.7062326     0.19956426   2.48   2.65   2.91
##      p_50    p_75     p_95     p_99   skewness  kurtosis     iqr
## 1    0.70    1.04     1.70     2.18  1.1173816  4.260117    0.64
## 2   61.80   62.50    63.80    65.60 -1.4218023 64.318726    1.50
## 3   57.00   59.00    61.00    64.00  0.7968737  5.801486    3.00
## 4 2401.00 5324.25 13107.10 17378.22  1.6183405  5.177360 4374.25
## 5    5.70    6.54     7.66     8.36  0.3786481  2.381801    1.83
## 6    5.71    6.54     7.65     8.34  2.4290393 94.079428    1.82
## 7    3.53    4.04     4.73     5.15  1.5325981 50.059515    1.13
##          range_98     range_80
## 1    [0.24, 2.18] [0.31, 1.51]
## 2    [57.9, 65.6]   [60, 63.3]
## 3        [53, 64]     [55, 60]
## 4 [429, 17378.22]  [646, 9821]
## 5    [4.02, 8.36] [4.36, 7.31]
## 6    [4.04, 8.34]  [4.36, 7.3]
## 7    [2.48, 5.15] [2.69, 4.52]

df_status(diamonds)

##    variable q_zeros p_zeros q_na p_na q_inf p_inf    type unique
## 1     carat       0    0.00    0    0     0     0 numeric    274
## 2       cut       0    0.00    0    0     0     0  factor      9
## 3     color       0    0.00    0    0     0     0  factor     10
## 4   clarity       0    0.00    0    0     0     0  factor      8
## 5     depth       0    0.00    0    0     0     0 numeric    185
## 6     table       0    0.00    0    0     0     0 numeric    127
## 7     price       0    0.00    0    0     0     0 integer  11603
## 8         x       8    0.01    1    0     0     0 numeric    554
## 9         y       8    0.01    0    0     0     0 numeric    554
## 10        z      20    0.04    0    0     0     0 numeric    376

Los resultados de df_status nos dice que hay ocho 0s en x e y, y viente 0s en z, si nos remitimos al diccionario de datos, nos damos cuenta que estas variables son distancia, por lo que no pueden ser 0. En este caso lo que haremos sera convertir esos 0s en NAs:

diamonds <- diamonds %>% 
  mutate(x = na_if(x, 0)) %>%
  mutate(y = na_if(y, 0)) %>%
  mutate(z = na_if(z, 0))

Con esto queda terminado el análisis estadístico.

Ahora para el análisis gráfico, analizaremos los boxplots e histogramas de cada variable para ver si podemos intuir sobre la posible existencia de valores atípicos. Para esto haremos un gráfico de facetas para los histogramas utilizando la función melt de reshape2 para hacer que la base sea larga:

diamonds %>% 
  select(is.numeric)  %>% 
  mutate(id = row_number()) %>% 
  melt(id.vars = "id") %>%
  ggplot(aes(x=value)) + 
  geom_histogram(bins = 50, fill = "purple") + 
  facet_wrap(~variable, scale = "free_x") +
  theme_bw()

De la misma manera que hicimos ese histograma, haremos el boxplot:

diamonds %>% 
  select(is.numeric)  %>% 
  mutate(id = row_number()) %>% 
  melt(id.vars = "id") %>%
  ggplot(aes(x=id, y=value)) +
  geom_boxplot(col = "black", fill = "purple", alpha = 0.7) +
  facet_wrap(~variable, scale = "free_y")

Con estos gráficos podemos inferir entonces que hay valores atípicos, los que marca el final del análisis gráfico y comenzamos a ver los atípicos.

Para el análisis de valores atípicos utilizaremos los umbrales de tukey que es una formula toda maluca, pero afortunadamente, R hace eso. funModeling, nos dirá cuales son las observaciones que están fuera de esta umbral.

En el boxplot y en el histograma se ve que la variable price tiene resto de atípicos, entonces usemos esta función para verlos:

tukey_outlier(diamonds$price)

## bottom_threshold    top_threshold 
##        -12172.75         18447.00

Esto nos da los umbrales del precio y ahora para ver cuales son y cuantas son las observaciones fuera de este umbral:

which(diamonds$price > tukey_outlier(diamonds$price)[2])

##   [1] 27601 27602 27603 27604 27605 27606 27607 27608 27609 27610 27611 27612
##  [13] 27613 27614 27615 27616 27617 27618 27619 27620 27621 27622 27623 27624
##  [25] 27625 27626 27627 27628 27629 27630 27631 27632 27633 27634 27635 27636
##  [37] 27637 27638 27639 27640 27641 27642 27643 27644 27645 27646 27647 27648
##  [49] 27649 27650 27651 27652 27653 27654 27655 27656 27657 27658 27659 27660
##  [61] 27661 27662 27663 27664 27665 27666 27667 27668 27669 27670 27671 27672
##  [73] 27673 27674 27675 27676 27677 27678 27679 27680 27681 27682 27683 27684
##  [85] 27685 27686 27687 27688 27689 27690 27721 27722 27723 27724 27725 27726
##  [97] 27727 27728 27729 27730 27731 27732 27733 27734 27735 27736 27737 27738
## [109] 27739 27740 27741 27742 27743 27744 27745 27746 27747 27748 27749 27750

length(which(diamonds$price > tukey_outlier(diamonds$price)[2]))

## [1] 120

Podemos ver entonces que hay 120 valores atípicos, es importantes entonces ver si fue un error de digitación o en el registro de información, si no es ninguno es mejor dejar esos errores en la muestra. Se puede hacer lo mismo para las otras variables pero es lo mismo.

Con esto hemos terminado las variables cuantitativas, por lo que empezaremos el análisis de las variables cualitativas, en el que se realizan los siguientes pasos:

Tablas de Frecuencia
Análisis Gráfico

Las tablas de frecuencia se hacen exactamente igual a como se hacen con las cuantitativas, es decir, con la función de describe:

describe(diamonds[v_cuali])

## diamonds[v_cuali] 
## 
##  3  Variables      53940  Observations
## --------------------------------------------------------------------------------
## cut 
##        n  missing distinct 
##    53940        0        9 
##                                                                       
## Value           Fair      good      Good     Ideal   premium   Premium
## Frequency       1610         1      4905     21551         1     13788
## Proportion     0.030     0.000     0.091     0.400     0.000     0.256
##                                         
## Value       Premiumm  Premmium Very Good
## Frequency          1         1     12082
## Proportion     0.000     0.000     0.224
## --------------------------------------------------------------------------------
## color 
##        n  missing distinct 
##    53940        0       10 
##                                                                       
## Value       best     D     E     F     g     G     H     I     J worst
## Frequency      1  6774  9797  9542     1 11291  8304  5422  2807     1
## Proportion 0.000 0.126 0.182 0.177 0.000 0.209 0.154 0.101 0.052 0.000
## --------------------------------------------------------------------------------
## clarity 
##        n  missing distinct 
##    53940        0        8 
##                                                           
## Value         I1    IF   SI1   SI2   VS1   VS2  VVS1  VVS2
## Frequency    741  1790 13065  9194  8171 12258  3655  5066
## Proportion 0.014 0.033 0.242 0.170 0.151 0.227 0.068 0.094
## --------------------------------------------------------------------------------

Podemos de primera vista que la palabra premium se escribio mal varias veces, por lo que es necesario corregir eso:

diamonds <- diamonds %>%
  mutate(cut = as.factor(str_to_sentence(cut))) #poner mayusculas en la primera letra

diamonds$cut <- diamonds$cut %>%
  recode(Premiumm = "Premium", Premmium = "Premium") #renombrar observaciones

De igual manera, la variable cut es una variable ordinal por lo que debemos hacer eso:

diamonds$cut <- ordered(diamonds$cut, levels = c("Fair", "Good", "Very good", 
                                                 "Premium", "Ideal"))

Con esto terminamos con la variable cut, pero la variable color también tiene un problema por lo que vamos a resolverlo:

diamonds$color <- diamonds$color %>%
  recode(best = "J", worst = "D")

diamonds <- diamonds %>%
  mutate(color = str_trim(color), #Quita los espacios que diferencian las observaciones
         color = str_to_upper(color)) #Pone todo el texto en mayuscula

Con esto, terminamos las tablas de frecuencia y de arreglar las observaciones. Ahora con eso terminamos empezamos el análisis gráfico que consta de realizar un gráfico de barras para las variables:

diamonds %>% 
  select(-is.numeric) %>% 
  mutate(id = row_number()) %>% 
  melt(id.vars = 'id') %>% 
  ggplot(aes(x = value)) + 
  geom_bar(fill = "purple", col = "black") + 
  facet_wrap(~ variable, scale = "free_x")

Con eso ya terminamos con los análisis de las variables, y llegamos a la etapa final, que es el análisis de la base en general, esto consiste de solo dos cosas:

Datos perdidos
Duplicados

Los datos perdidos son los NAs en la base, entonces para crear una base aparte sin los NAs debemos primero identificarlos:

table(is.na(diamonds))

## 
##  FALSE   TRUE 
## 539363     37

colSums(is.na(diamonds))

##   carat     cut   color clarity   depth   table   price       x       y       z 
##       0       0       0       0       0       0       0       9       8      20

diamonds[!complete.cases(diamonds),]

##       carat       cut color clarity depth table price    x    y    z
## 72     0.30 Very good     H     SI1  63.1    56   554 4.29   NA 2.70
## 73     0.30   Premium     H     SI1  62.9    59   554   NA 4.24 2.68
## 2208   1.00   Premium     G     SI2  59.1    59  3142 6.55 6.48   NA
## 2315   1.01   Premium     H      I1  58.1    59  3167 6.66 6.60   NA
## 4792   1.10   Premium     G     SI2  63.0    59  3696 6.50 6.47   NA
## 5472   1.01   Premium     F     SI2  59.2    58  3837 6.50 6.47   NA
## 10168  1.50      Good     G      I1  64.0    61  4731 7.15 7.04   NA
## 11183  1.07     Ideal     F     SI2  61.6    56  4954   NA 6.62   NA
## 11964  1.00 Very good     H     VS2  63.3    53  5139   NA   NA   NA
## 13602  1.15     Ideal     G     VS2  59.2    56  5564 6.88 6.83   NA
## 15952  1.14      Fair     G     VS1  57.5    67  6381   NA   NA   NA
## 24395  2.18   Premium     H     SI2  59.4    61 12631 8.49 8.45   NA
## 24521  1.56     Ideal     G     VS2  62.2    54 12800   NA   NA   NA
## 26124  2.25   Premium     I     SI1  61.3    58 15397 8.52 8.42   NA
## 26244  1.20   Premium     D    VVS1  62.1    59 15686   NA   NA   NA
## 27113  2.20   Premium     H     SI1  61.2    59 17265 8.42 8.37   NA
## 27430  2.25   Premium     H     SI2  62.8    59 18034   NA   NA   NA
## 27504  2.02   Premium     H     VS2  62.7    53 18207 8.02 7.95   NA
## 27740  2.80      Good     G     SI2  63.8    58 18788 8.90 8.85   NA
## 49557  0.71      Good     F     SI2  64.1    60  2130   NA   NA   NA
## 49558  0.71      Good     F     SI2  64.1    60  2130   NA   NA   NA
## 51507  1.12   Premium     G      I1  60.4    59  2383 6.71 6.67   NA

De igual manera podemos hacer un plot para visualizar estos valores perdidos con la librería de finalfit, sin embargo, con esta base de diamonds no hay suficientes NAs para visualizarlos por lo que usaremos otra para mostrar el ejemplo:

missing_plot(colon_s) #base de finalfit

Ahora si podemos eliminar los NAs en una base aparte:

diamonds_sinna<- na.omit(diamonds)

Como ultimo paso, identificaremos los duplicados en la base para eliminarlos, es importante reconocer que existen diferentes tipos de duplicados, totales que son duplicados en todas las variables, o parciales que hay una variables que no están duplicadas. Entonces miramos los duplicados en la base:

sum(duplicated(diamonds))

## [1] 146

head(diamonds[duplicated(diamonds),])

##      carat   cut color clarity depth table price    x    y    z
## 1006  0.79 Ideal     G     SI1  62.3    57  2898 5.90 5.85 3.66
## 1007  0.79 Ideal     G     SI1  62.3    57  2898 5.90 5.85 3.66
## 1008  0.79 Ideal     G     SI1  62.3    57  2898 5.90 5.85 3.66
## 1009  0.79 Ideal     G     SI1  62.3    57  2898 5.90 5.85 3.66
## 2026  1.52  Good     E      I1  57.3    58  3105 7.53 7.42 4.28
## 2184  1.00  Fair     E     SI2  67.0    53  3136 6.19 6.13 4.13

Con los duplicados ya identificados, los podemos eliminar, ya que es preferible que las bases no tengan duplicados:

diamonds_sindup <- distinct(diamonds)

Con esto finalizando, volvemos a limpiar el workspace.

Unidad 5: Clustering

La tarea de clustering consiste en la agrupación de individuos según características similares, para esto se toman en cuenta dos cosas, maximizar la distancia inter-cluster y minimizar la distancia intra-cluster, es decir, que los grupos sean diferente entre ellos, pero dentro sean similares. El proceso de clustering tiene los pasos para determinar los grupos:

Medidas de Similitud: se utilizan para determinar la cercania de los individuos, hay métricas de distancia, de correlación y de asociación. Nosotros nos centraremos en distancia, específicamente en la euclidiana
Algoritmos para la Formación de Clústeres: ya después de determinar la distancia de los individuos, se debe decidir que algoritmo se utilizara para formar los clústeres. Existen dos tipos: jerárquicos y particionados. Los jerárquicos construyen subconjuntos de clústeres organizados en jerarquicamente en un árbol (dendrograma) y los particionados forman grupos que son mutuamente excluyentes.
Criterio para Determinar el Número de Clústeres: consiste en elegir el numero de grupos razonables, para esto existen bastantes metodos que veremos cuando entremos en materia.

Antes de iniciar el clave aclarar que existen dos tipos de clústeres jerárquicos, los aglomerativos, que comienza con cada observación siendo su propio clúster y se va expandiendo, y los de división que son lo opuesto.

Ahora si, vamos a cargar los paquetes y bases para esta unidad:

require(pacman)
p_load(tidyverse, NbClust, factoextra, cluster, rio)

datos_originales <- import("bases/Credit_Card_Customer_Data.csv")

Como dijimos antes, primero debemos ver los datos y asegurarnos que no hallan problemas:

glimpse(datos_originales)

## Rows: 660
## Columns: 7
## $ Sl_No               <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
## $ `Customer Key`      <int> 87073, 38414, 17341, 40496, 47437, 58634, 48370, 3…
## $ Avg_Credit_Limit    <int> 100000, 50000, 50000, 30000, 100000, 20000, 100000…
## $ Total_Credit_Cards  <int> 2, 3, 7, 5, 6, 3, 5, 3, 2, 4, 4, 3, 1, 1, 2, 2, 2,…
## $ Total_visits_bank   <int> 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 2, 0, 1,…
## $ Total_visits_online <int> 1, 10, 3, 1, 12, 1, 11, 1, 2, 1, 5, 2, 5, 5, 4, 5,…
## $ Total_calls_made    <int> 0, 9, 4, 4, 3, 8, 2, 1, 2, 7, 5, 7, 9, 6, 6, 7, 4,…

Vemos que las primeras dos variables no nos sirven puesto que son de identificación, por lo que las eliminaremos:

datos <- datos_originales[,-c(1,2)]
glimpse(datos)

## Rows: 660
## Columns: 5
## $ Avg_Credit_Limit    <int> 100000, 50000, 50000, 30000, 100000, 20000, 100000…
## $ Total_Credit_Cards  <int> 2, 3, 7, 5, 6, 3, 5, 3, 2, 4, 4, 3, 1, 1, 2, 2, 2,…
## $ Total_visits_bank   <int> 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 2, 0, 1,…
## $ Total_visits_online <int> 1, 10, 3, 1, 12, 1, 11, 1, 2, 1, 5, 2, 5, 5, 4, 5,…
## $ Total_calls_made    <int> 0, 9, 4, 4, 3, 8, 2, 1, 2, 7, 5, 7, 9, 6, 6, 7, 4,…

Ya con eso listo, podemos empezar estandarizando los datos, esto significa dejar los datos en una escala similar, puesto que como estamos hablando de una distancia, debemos tener los datos de manera similar, para esto usamos scale:

datos_est <- as.data.frame(scale(datos))
glimpse(datos_est)

## Rows: 660
## Columns: 5
## $ Avg_Credit_Limit    <dbl> 1.7388680, 0.4099816, 0.4099816, -0.1215730, 1.738…
## $ Total_Credit_Cards  <dbl> -1.2482780, -0.7869883, 1.0581707, 0.1355912, 0.59…
## $ Total_visits_bank   <dbl> -0.8597985, -1.4726139, -0.8597985, -0.8597985, -1…
## $ Total_visits_online <dbl> -0.5470748, 2.5186084, 0.1341882, -0.5470748, 3.19…
## $ Total_calls_made    <dbl> -1.2505889, 1.8904250, 0.1454173, 0.1454173, -0.20…

Ahora, debemos calcular la matriz de distancia, que es lo que nos permitirá saber la distancia entre los puntos, para construir el dendrograma y poder realizar los clústeres:

datos_dist <- dist(datos_est)

Ahora, para realizar el dendrograma utilizamos la función de hclust con el método de centroide, que es uno de los varios métodos de hacer un clustering jerárquico:

HCA_centroide <- hclust(datos_dist, method = "centroid")
HCA_centroide

## 
## Call:
## hclust(d = datos_dist, method = "centroid")
## 
## Cluster method   : centroid 
## Distance         : euclidean 
## Number of objects: 660

Y para ver el dendrograma generado realizamos la siguiente función, antes es necesario recalcar que como son varias observaciones, no se va a ver bien:

plot(HCA_centroide, main = "Dendrograma método centroide", hang = -1, 
     cex = 0.2, ylab = "Distancia")

Con eso definido, debemos ahora escoger el numero optimo de clústeres, para esto emplearemos la función NbClust del paquete que lleva el mismo nombre, esto nos ayudara a 30 indices para la decisión:

res_centroid <- NbClust(datos_est, distance = "euclidean", min.nc = 2, 
                        max.nc = 10, method = "centroid", index = "all")

## [1] "Frey index : No clustering structure in this data set"

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
##

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 6 proposed 2 as the best number of clusters 
## * 2 proposed 3 as the best number of clusters 
## * 1 proposed 4 as the best number of clusters 
## * 12 proposed 5 as the best number of clusters 
## * 2 proposed 10 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  5 
##  
##  
## *******************************************************************

Para visualizar los números de clústeres que cada indice recomendó realizamos:

votes_centroid <- as.data.frame(table(res_centroid$Best.nc[1, ])) %>%
  rename(clusters = Var1, votes = Freq) %>%
  mutate(clusters = as.numeric(as.character(clusters)))

votes_centroid %>%
  ggplot(aes(x = clusters, y = votes)) +
  geom_col(fill = "purple", col = "black") +
  geom_label(aes(label = votes)) +
  labs(title = "Número óptimo de clusters según NbClust",
       x = "Número de clusters",
       y = "Cantidad de índices que lo recomiendan") +
  theme_minimal()

Ahora vamos a utilizar un criterio más para la decisión del clúster optimo y este es el uso de la función silhouette en conjunto con cutree para visualizar las siluetas de los clústeres y ver que tan bien posicionados están los individuos. Para interpretar esta silueta debemos tener en cuenta estos rangos:

No hay estructura clara: menor o igual a 0.25
Estructura débil: entre 0.26 y 0.50
Razonablemente bueno: entre 0.51 y 0.70
Muy fuerte: mayor a 0.70

Con eso vamos a evaluar las siluetas de 5, 4, 3 y 2:

plot(silhouette(cutree(HCA_centroide, 5), datos_dist), border = NA)

plot(silhouette(cutree(HCA_centroide, 4), datos_dist), border = NA)

plot(silhouette(cutree(HCA_centroide, 3), datos_dist), border = NA)

plot(silhouette(cutree(HCA_centroide, 2), datos_dist), border = NA)

Viendo estos resultados, podemos concluir que el numero de clústeres optimo es 2 si utilizamos el método centroide:

grupos_HCA_aglo <- cutree(HCA_centroide ,2)

Ahora visualicemos los clústeres y el dendrograma. Primero el dendrograma y luego los clústeres:

plot(HCA_centroide, main = "Dendrograma método del centroide", 
     hang = -1, cex = 0.1, ylab = "Distancia")
rect.hclust(HCA_centroide, k = 2, border = 2:5 )

fviz_cluster(list(data = datos_est, cluster = grupos_HCA_aglo))+
  theme_bw() + 
  scale_colour_manual(values = c("#2700ff","#a200ff")) +
  scale_fill_manual(values = c("#2700ff","#a200ff"))

Ya con eso vimos que son los clústeres jerárquicos, pero mencione que existían también los particionados. Para este tipo de clústeres utilizamos el aprendizaje de maquina no supervisado para realizar los clústeres, esto significa que debemos utilizar semillas para que el aspecto aleatorio del algoritmo no cambie el resultado cada que sea utilizado. Para empezar como ya tenemos los datos estandarizados saltamos directo al análisis de los clústeres óptimos que se realiza de la misma manera:

set.seed(171124)

res_kmeans <- NbClust(datos_est, distance = "euclidean",
                      min.nc = 2, max.nc = 10, 
                      method = "kmeans", index = "all")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
##

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 1 proposed 2 as the best number of clusters 
## * 20 proposed 3 as the best number of clusters 
## * 1 proposed 6 as the best number of clusters 
## * 1 proposed 10 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************

votes_kmeans <- as.data.frame(table(res_kmeans$Best.nc[1, ])) %>%
  rename(clusters = Var1, votes = Freq) %>%
  mutate(clusters = as.numeric(as.character(clusters)))

votes_kmeans %>%
  ggplot(aes(x = clusters, y = votes)) +
  geom_col(fill = "purple", col = "black") +
  geom_label(aes(label = votes)) +
  labs(title = "Número óptimo de clusters según NbClust",
       x = "Número de clusters",
       y = "Cantidad de índices que lo recomiendan") +
  theme_minimal()

set.seed(171124)
kmeans3 = kmeans(datos_est, 3)
plot(silhouette(kmeans3$cluster, dist(datos_est)), border = NA)

set.seed(171124)
kmeans2 = kmeans(datos_est, 2)
plot(silhouette(kmeans2$cluster, dist(datos_est)), border = NA)

Una vez más, utilizando un método de clustering particionado la cantidad de clústeres óptimos es 2:

cluster_kmeans <- kmeans(x = datos_est, centers = 2)

fviz_cluster(cluster_kmeans, datos_est)+
  theme_bw() + 
  scale_colour_manual(values = c("#2700ff","#a200ff")) +
  scale_fill_manual(values = c("#2700ff","#a200ff"))

Antes de finalizar esta unidad, se puede hacer algo extra que es realizar un análisis ANOVA y grafico de los clústeres, para esto proveeré las fotos y el código de un taller que realice donde utilice esto para analizar los clústeres:

Análisis Gráfico:

p1 <- taller5.clust %>%
  ggplot(aes(x = Alcohol, y = Cluster, col = Cluster, fill = Cluster)) +
  geom_boxplot(alpha = 0.5) +
  scale_color_manual(aesthetics = c("fill", "colour"), 
                     values = c("1" = "#a200ff", "2" = "#f700ff", 
                                "3" = "#2700ff", "4" = "#ff006c")) +
  labs(x = "Porcentaje de Alcohol", y = "Grupos") +
  theme_bw() +
  theme(legend.position = "none")

p2 <- taller5.clust %>%
  ggplot(aes(x = Malic_Acid, y = Cluster, col = Cluster, fill = Cluster)) +
  geom_boxplot(alpha = 0.5) +
  scale_color_manual(aesthetics = c("fill", "colour"), 
                     values = c("1" = "#a200ff", "2" = "#f700ff", 
                                "3" = "#2700ff", "4" = "#ff006c")) +
  labs(x = "Cantidad de Acido Malico", y = "Grupos") +
  theme_bw() +
  theme(legend.position = "none")

p3 <- taller5.clust %>%
  ggplot(aes(x = Ash, y = Cluster, col = Cluster, fill = Cluster)) +
  geom_boxplot(alpha = 0.5) +
  scale_color_manual(aesthetics = c("fill", "colour"), 
                     values = c("1" = "#a200ff", "2" = "#f700ff", 
                                "3" = "#2700ff", "4" = "#ff006c")) +
  labs(x = "Cantidad de Minerales", y = "Grupos") +
  theme_bw() +
  theme(legend.position = "none")

grid.arrange(p1, p2, p3, nrow = 3)

Tabla ANOVA:

tema <- ttheme_default(core = list(bg_params = list(fill = "white", 
                                                    col = "black")),
  colhead = list(bg_params = list(fill = "#a200ff", col = "black"),
    fg_params = list(col = "black", fontface = "bold")))

tabla1 <- aggregate(taller5.clust[1:4], list(taller5.clust$Cluster), mean)
names(tabla1) <- c("Grupo", "Alcohol", "Acido Malico", "Minerales", 
                   "Alcanilidad")
rownames(tabla1) <- NULL
tabla1.plot <- tableGrob(tabla1, rows = NULL, theme = tema)

tabla2 <- aggregate(taller5.clust[5:8], list(taller5.clust$Cluster), mean)
names(tabla2) <- c("Grupo", "Magnesio", "Fenólicos", "Flavonoides",
                      "No Flavonoides")
rownames(tabla2) <- NULL
tabla2.plot <- tableGrob(tabla2, rows = NULL, theme = tema)

tabla3 <- aggregate(taller5.clust[9:13], list(taller5.clust$Cluster), mean)
names(tabla3) <- c("Grupo", "Proantocianina", "Intensidad", "Tonalidad",
                      "Absorbancia", "Prolina")
rownames(tabla3) <- NULL
tabla3.plot <- tableGrob(tabla3, rows = NULL, theme = tema)

grid.arrange(tabla1.plot, tabla2.plot, tabla3.plot, nrow = 3)

Test del ANOVA: Una vez se realice la tabla, es necesario evaluar si las variables pueden entrar al análisis, para esto, se realizan las siguientes pruebas:
- ANOVA: se utiliza para verificar que las medias entre los grupos sean iguales, donde H0 es que son iguales, y que HA es que al menos una es diferente
- Levene: se utiliza para verificar la homocedasticidad (igualdad de varianzas), donde H0 es que son iguales y HA es que al menos una es diferente.
- Shapiro: se utiliza para verificar que los datos se distribuyan de manera normal, donde H0 es que se distribuyen de manera normal y HA es que no
- Oneway (ANOVA de Welch): se utiliza cuando no se cumple el supuesto de homocedasticidad, donde H0 es que las medias son iguales y HA es que al menos una media es significativamente diferente al resto, lo que hace que el análisis sea importante
Para estos test, recordemos que rechazamos H0 si el p-value es mejor que 0.01, y solo se realiza el ANOVA de Welch si rechazamos H0 en Levene. Rechazamos una variable por completo si no son normales y no hay homocedasticidad. En caso de que se cumpla la homocedasticidad pero no el principio de normalidad vemos el tamaño de la muestra, si esta es suficientemente grande:
- No se acepta: menor a 30
- Se puede relajar el principio: mayor o igual a 30
- Se puede ignorar con confianza en la mayoría de casos: entre 50 y 100
- Se puede rechazar por completo: mayor a 200
```
#ANOVA
aov.noflavonoides <- aov(Nonflavanoid_Phenols ~ Cluster, data = taller5.clust)
summary(aov.noflavonoides)

#Levene
leveneTest(aov.noflavonoides)

#Shapiro
residual.noflavonoides <- residuals(aov.noflavonoides)
shapiro.test(x = residual.noflavonoides)

#ANOVA de Welch
oneway.test(Nonflavanoid_Phenols ~ Cluster, data = taller5.clust)
```

Ya finalmente se pueden eliminar las variables que no cumplen con los principios de normalidad, de homocedasticidad, para esto se realiza la tabla de la misma manera, solo que no se incluyen aquellas variables que no sirven.

Una vez más eliminamos el workspace.

Unidad 6: Resumen

Otra tarea del BA es la tarea de resumir datos. En este caso nos vamos a enfocar en resumir textos por medio del uso de nubes de palabra, esto es una manera de analizar textos de manera gráfica utilizando la minería de textos.

Los paquetes y bases que utilizaremos para esta unidad serán:

require(pacman)
p_load(tidyverse, tm, wordcloud, wordcloud2, RColorBrewer)

discurso <- readLines("bases/Santos.txt",encoding = "UTF-8")

Con el discurso ya cargado, debemos usar herramientas de minería de texto que fue mencionado previamente. Primero, para poder realizar esto debemos crear un objeto que pueda ser manipulado de esta manera, para esto utilizaremos la librería tm:

texto.Santos <- Corpus(VectorSource(discurso))

Ahora debemos realizar las siguientes tareas:

Pasar todas las palabras a minúsculas
Eliminar los espacios en blanco
Quitar los signos de puntuación
Borrar los numeros
Remover palabras vacías como los, las, tendremos, etc.

Entonces haremos esto en ese mismo orden:

texto.Santos <- tm_map(texto.Santos, tolower)

texto.Santos <- tm_map(texto.Santos, stripWhitespace)

texto.Santos <- tm_map(texto.Santos, removePunctuation)

texto.Santos <- tm_map(texto.Santos, removeNumbers)

texto.Santos <- tm_map(texto.Santos, removeWords, stopwords("spanish"))

Ahora con el texto limpio, debemos convertirlo en una base de datos para poder calcular las frecuencias de aparición de cada palabra:

Santos <- TermDocumentMatrix(texto.Santos)

Antes de pasar a la base de datos podemos realizar unas cosas como encontrar las palabras de alta y baja frecuencia, y correlaciones entre palabras, de esta manera:

findFreqTerms(Santos, lowfreq = 10)

## [1] "paz"      "acuerdo"  "hoy"      "colombia" "guerra"   "gracias"

findAssocs(Santos, "colombia", 0.5)

## $colombia
##   educada   equidad   felices   permita progresar       ser 
##      0.58      0.58      0.58      0.58      0.58      0.58

Se pueden hallar la asociación de varias y diferentes indices de asociación simplemente haciendo una lista c().

Ahora si ya podemos crear la base de datos:

Santos.matriz <- as.matrix(Santos)

data.Santos <- sort(rowSums(Santos.matriz), decreasing = TRUE)
data.Santos <- data.frame(word = names(data.Santos),freq = data.Santos)

Ya con esto finalizado podemos visualizar el texto, esto se puede realizar de tres maneras, utilizando un gráfico de barras, de puntos o la nube:

ggplot(data.Santos[1:10, ], aes(x = reorder(word, -freq), y = freq)) +
  geom_col(col = "black", fill = "purple") +
  geom_label(aes(label = freq)) +
  coord_flip() +
  labs(title = "10 Palabras con Mayor Frecuencia en el Discurso de Santos",
       y = "Frecuencia", x = "Palabra", caption = "Fuente: Elaboración Propia") +
  theme_bw()

ggplot(data.Santos[1:10, ], aes(x = reorder(word, -freq), y = freq)) +
  geom_point(col = "purple", size = 3) +
  geom_label(aes(label = freq), hjust = 0) +
  coord_flip() +
  labs(title = "10 Palabras con Mayor Frecuencia en el Discurso de Santos",
       y = "Frecuencia", x = "Palabra", caption = "Fuente: Elaboración Propia") +
  theme_bw()

wordcloud(words = data.Santos$word, freq = data.Santos$freq, min.freq = 1,
          random.order = FALSE, colors = brewer.pal(10, "Blues"))

wordcloud2(data = data.Santos, size = 1.1, color = brewer.pal(9, "Blues"),
           maxRotation = pi/2, shape = 'circle', ellipticity = 0.2)

Con eso hecho, es posible comparar dos textos utilizando estas nubes de palabras, de esta manera. Primero creamos y manipulamos el discurso de timochenko:

discurso1 <- readLines("bases/Timochenko.txt",encoding = "UTF-8")

texto.Timochenko <- Corpus(VectorSource(discurso1))

texto.Timochenko <- tm_map(texto.Timochenko, tolower)
texto.Timochenko <- tm_map(texto.Timochenko, stripWhitespace)
texto.Timochenko <- tm_map(texto.Timochenko, removePunctuation)
texto.Timochenko <- tm_map(texto.Timochenko, removeNumbers)
texto.Timochenko <- tm_map(texto.Timochenko, removeWords, stopwords("spanish"))

Timochenko <- TermDocumentMatrix(texto.Timochenko)

Timochenko.matriz <- as.matrix(Timochenko)

data.Timochenko <- sort(rowSums(Timochenko.matriz), decreasing = TRUE)
data.Timochenko <- data.frame(word = names(data.Timochenko),freq = data.Timochenko)

Para realizar la comparación debemos entonces crear una base que tenga ambos discursos, organizarla y crear la nube de comparación:

names(data.Santos)[2] <- "Santos"
names(data.Timochenko)[2] <- "Timochenko"

data.completa <- full_join(data.Santos, data.Timochenko)

## Joining with `by = join_by(word)`

#Existen palabras que hay en discurso que el otro no tiene
sum(is.na(data.completa$Santos))

## [1] 853

sum(is.na(data.completa$Timochenko))

## [1] 413

data.completa$Santos <- data.completa$Santos %>% 
  replace_na(0)
data.completa$Timochenko <- data.completa$Timochenko %>% 
  replace_na(0)

row.names(data.completa) = data.completa[,1]
data.completa <- as.matrix(data.completa[,-1])

comparison.cloud(data.completa)

De igual manera se puede realizar nubes de palabras en forma de letra o palabra utilizando la función letterCloud o utilizar formas unicas usando el argumento de figPath en la función wordcloud2 donde se introduce la dirección del png, de la siguiente manera:

#Letra/Palabra
names(data.Santos)[2] <- "freq"
letterCloud(data.Santos,
size = 0.9, word = "PAZ", color= "skyblue")

#Imagen
wordcloud2(data.Santos, figPath = "peaceAndLove.jpg", size = 1.5,
color = "skyblue", backgroundColor = "black", shuffle = FALSE)

Borramos workspace y seguimos.

Tercer Corte

Unidad 7: Análisis de Redes

El análisis de redes surge de la necesidad de entender la interdependencia y las conexiones informales dentro de la sociedad, los grafos son gráficos que conectan a personas y son utilizados precisamente para analizar este tipo de relaciones.

Estas conexiones están definidas por una matriz de adyacente cuadrada donde la diagonal debe ser 0 y el numero en la matriz representa la cantidad de conexiones.

Antes de empezar a entender el procedimiento para realizar un grafo, debemos aprender los dos tipos de grafos:

Ponderado: el grosor de la conexión depende de la cantidad de conexiones que la persona tiene con el otro
Dirigido: La conexión tiene una flecha que muestra hacia que lado va la conexión

El grafo puede ser ponderado y dirigido, ponderado y no dirigido, no ponderado y no dirigido o no ponderado y no dirigido.

Ahora si para poder realizar los grafos necesitaremos los siguientes paquetes y la siguiente base:

require(pacman)
p_load(tidyverse, igraph, reshape2, igraph)

Croacia <- import("bases/partidos.xlsx", sheet = "Croacia")
head(Croacia)

##   Subasic Vrsaljko Strinic Perisic Lovren Rakitic Modirc Brozovic Mandzukic
## 1      NA        1      NA      NA      7       1      1        1         1
## 2       1       NA      NA       5     15       4     10       10         5
## 3      NA       NA      NA       6     NA       5     NA        4         1
## 4      NA        2       3      NA     NA       3      8       NA         1
## 5       3       10      NA       1     NA       8      9       17         3
## 6       1        2       6       5      5      NA      6        9         2
##   Rebic Vida Kramaric Pjaca
## 1    NA    3       NA    NA
## 2     2   NA       NA     1
## 3     1    1       NA    NA
## 4     1    2       NA    NA
## 5    NA    7       NA    NA
## 6     7    5        1    NA

Lo primero como siempre es limpiar los datos, y en este caso vemos hay muchos NAs, por lo primero que realizaremos sera convertir estos NAs en 0s y convertir la base en una matriz:

Croacia[is.na(Croacia)] <- 0
Croacia <- as.matrix(Croacia)

Con esto ya podemos hacer el grafo. Realizaremos el grafo de las 4 maneras que se menciono previamente:

#No ponderado, no dirigido
Croacia.Und.graph <- graph.adjacency(Croacia, mode = "undirected", diag = FALSE)

V(Croacia.Und.graph)$color <- "green" #Color del objeto
V(Croacia.Und.graph)$shape <- "sphere" #Forma del objeto
E(Croacia.Und.graph)$color <- "gray" #Color de las lineas

plot(Croacia.Und.graph)

#No ponderado, dirigido
Croacia.Dir.graph <- graph.adjacency(Croacia, mode = "directed", diag = FALSE)

V(Croacia.Dir.graph)$color <- "green"
V(Croacia.Dir.graph)$shape <- "sphere"
E(Croacia.Dir.graph)$color <- "gray"

plot(Croacia.Dir.graph, edge.arrow.size = 0.8)

#Ponderado, no dirigido
Croacia.Und.Pond.graph <- graph.adjacency(Croacia, weighted = TRUE, 
                                          mode ="undirected", diag = FALSE)

V(Croacia.Und.Pond.graph)$color <- "green"
V(Croacia.Und.Pond.graph)$shape <- "sphere"
E(Croacia.Und.Pond.graph)$color <- "gray"

plot(Croacia.Und.Pond.graph, edge.width = E(Croacia.Und.Pond.graph)$weight)

#Ponderado, dirigido
Croacia.Dir.Pond.graph <- graph.adjacency(Croacia, weighted = TRUE, 
                                          mode ="directed", diag = FALSE)

V(Croacia.Dir.Pond.graph)$color <- "green"
V(Croacia.Dir.Pond.graph)$shape <- "sphere"
E(Croacia.Dir.Pond.graph)$color <- "gray"

plot(Croacia.Dir.Pond.graph, edge.width = E(Croacia.Dir.Pond.graph)$weight, 
     edge.arrow.size = 0.8)

Ahora bien, con el grafo no tenemos suficientes información sobre la red, incluso, como se puedo notar, es bastante complicado poder analizar eso por el hecho de que la lineas se sobreponen. Para esto utilizaremos indicadores de centralidad que nos darán más información sobre la red, estos vienen de varios sabores.

Primero están las medidas locales:

Grado de Centralidad: el numero de enlaces relacionados con cada nodo
Betweenness: el numero de veces que un nodo actúa como un puente entre dos nodos

degree(Croacia.Und.graph, mode = "all")

##   Subasic  Vrsaljko   Strinic   Perisic    Lovren   Rakitic    Modirc  Brozovic 
##        15        69        27        40        68        69        87       104 
## Mandzukic     Rebic      Vida  Kramaric     Pjaca 
##        35        23        44        17         4

betweenness(Croacia.Und.graph)

##     Subasic    Vrsaljko     Strinic     Perisic      Lovren     Rakitic 
## 0.000000000 2.336537714 0.074299223 0.408524883 0.175000000 2.278725916 
##      Modirc    Brozovic   Mandzukic       Rebic        Vida    Kramaric 
## 4.401537430 6.109251656 0.433124848 0.009433962 1.253863947 0.519700420 
##       Pjaca 
## 0.000000000

Luego están las medidas globales:

Closeness: calcula las rutas más cortas entre nodos, es decir, mide cuantos pasos requieren para conectarse uno de los nodos desde otro determinado.
Puntaje de Autoridad: mide la importancia relativa de cada nodo

closeness(Croacia.Und.graph, mode = "all")

##    Subasic   Vrsaljko    Strinic    Perisic     Lovren    Rakitic     Modirc 
## 0.05882353 0.07692308 0.06250000 0.07142857 0.06250000 0.07692308 0.08333333 
##   Brozovic  Mandzukic      Rebic       Vida   Kramaric      Pjaca 
## 0.08333333 0.07142857 0.06250000 0.07692308 0.06250000 0.05000000

authority.score(Croacia.Und.graph)$vector

##    Subasic   Vrsaljko    Strinic    Perisic     Lovren    Rakitic     Modirc 
## 0.15694159 0.77719921 0.25573956 0.37877595 0.73846071 0.64157987 0.88221620 
##   Brozovic  Mandzukic      Rebic       Vida   Kramaric      Pjaca 
## 1.00000000 0.36573244 0.22360943 0.44623143 0.16460428 0.04355986

También hay otras medidas que no entran a ninguna de las dos categorías previas:

Distancia: la cantidad de aristas en la ruta más corta que conecta un par de nodos
Excentricidad: la distancia de un nodo particular al más lejano

mean_distance(Croacia.Und.graph)

## [1] 1.230769

eccentricity(Croacia.Und.graph, mode = "all")

##   Subasic  Vrsaljko   Strinic   Perisic    Lovren   Rakitic    Modirc  Brozovic 
##         2         2         2         2         2         2         1         1 
## Mandzukic     Rebic      Vida  Kramaric     Pjaca 
##         2         2         2         2         2

Por ultimo, tenemos las medidas agregadas de la red:

Diámetro: mide la máxima excentricidad
Densidad: la proporción de pares de nodos conectados en la red sobre todos los posibles (porcentaje de conexiones efectivas)
Transitividad: mide la ocurrencia de pequeñas subredes de tres nodos totalmente conexas
Reciprocidad: mide la probabilidad de que si existe una conexión sale de un nodo x a uno y, haya otro que salga de y a x

diameter(Croacia.Und.graph)

## [1] 2

Croacia.1.0 <- Croacia
for (i in 1:ncol(Croacia.1.0)){
  for (j in 1:nrow(Croacia.1.0)){
    Croacia.1.0[j,i]<-ifelse(Croacia.1.0[j,i]>0,1,Croacia.1.0[j,i])
  }
}
GrafDens <- graph.adjacency(Croacia.1.0, mode = "undirected", weighted = NULL,
                            diag = FALSE, add.colnames=NULL)
edge_density(GrafDens)

## [1] 0.7692308

transitivity(Croacia.Und.graph)

## [1] 0.8269962

reciprocity(Croacia.Und.graph)

## [1] 1

Ya que acabamos de ver como interpretar un grafo, es útil saber que hay otras maneras de representar grafos, que es por medio del uso de heatmaps:

row.names(Croacia) <- variable.names(Croacia)
Croacia.por.parejas <- melt(Croacia)

ggplot(data = Croacia.por.parejas, aes(x = Var2, y = Var1, fill=value)) + 
  geom_tile() +
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
  scale_fill_gradient(low = "white", high = "purple") +
  xlab("receptor") + ylab("pasador") + labs(fill = "pases")

Y por ultimo también es importante mencionar que es posible encontrar clústers en las redes:

HC.Croacia.Und.Pond.graph <- cluster_edge_betweenness(Croacia.Und.Pond.graph)
#Halla los clusters

dendPlot(HC.Croacia.Und.Pond.graph, mode = "hclust", rect = 2) #Visualiza el dendrograma

#Visualizamos los clusters en el grafo
plot(HC.Croacia.Und.Pond.graph, Croacia.Und.Pond.graph)

length(HC.Croacia.Und.Pond.graph) #Cantidad de clusters

## [1] 2

membership(HC.Croacia.Und.Pond.graph)

##   Subasic  Vrsaljko   Strinic   Perisic    Lovren   Rakitic    Modirc  Brozovic 
##         1         2         2         2         1         2         2         2 
## Mandzukic     Rebic      Vida  Kramaric     Pjaca 
##         2         2         1         2         2

Limpiamos el workspace.

Unidad 8: Asociaciones (MBA)

El análisis de canasta de mercado es un conjunto de técnicas de minería de datos que usan los minoristas para encontrar asociaciones en sus clientes para mejorar sus ventas, es decir, es la tarea de asociación. Este análisis es un tipo de aprendizaje de maquina no supervisado donde se usan algoritmos y datos de transacciones, también conocida como tirillas de compra. Para este análisis existen dos partes:

Antecedente: el producto que causa la compra del otro
Consecuente: el producto resultado de la asociación

Ahora si entrando en materia de código. Cargaremos los paquetes necesarios y la base de datos:

require(pacman)
p_load(tidyverse, arules, arulesViz, rio)

datos <- import("bases/Groceries_dataset2025.csv")
glimpse(datos)

## Rows: 38,600
## Columns: 3
## $ Member_number   <int> 1808, 2552, 2300, 1187, 3037, 4941, 4501, 3803, 2762, …
## $ Date            <chr> "21-07-2015", "05-01-2015", "19-09-2015", "12-12-2015"…
## $ itemDescription <chr> "tropical fruit", "whole milk", "pip fruit", "other ve…

Antes de empezar, es necesario que los datos estén en formato de transacción, lo que implica que tendremos que agrupar los ítems por el ID de la transacción, esto nos dará los productos comprados en una canasta. Antes de poder convertirlos en una transacción tienen que ser una lista, para esto usamos la función split() de la base de R:

datos.lista <- split(datos$itemDescription, datos$Member_number)
class(datos.lista)

## [1] "list"

head(datos.lista)

## $`1000`
##  [1] "soda"                "canned beer"         "sausage"            
##  [4] "sausage"             "whole milk"          "whole milk"         
##  [7] "pickled vegetables"  "misc. beverages"     "semi-finished bread"
## [10] "hygiene articles"    "yogurt"              "pastry"             
## [13] "salty snack"        
## 
## $`1001`
##  [1] "frankfurter"        "frankfurter"        "beef"              
##  [4] "sausage"            "whole milk"         "soda"              
##  [7] "curd"               "white bread"        "whole milk"        
## [10] "soda"               "whipped/sour cream" "rolls/buns"        
## 
## $`1002`
## [1] "tropical fruit"      "butter milk"         "butter"             
## [4] "frozen vegetables"   "sugar"               "specialty chocolate"
## [7] "whole milk"          "other vegetables"   
## 
## $`1003`
## [1] "sausage"         "root vegetables" "rolls/buns"      "detergent"      
## [5] "frozen meals"    "rolls/buns"      "dental care"     "rolls/buns"     
## 
## $`1004`
##  [1] "other vegetables"          "pip fruit"                
##  [3] "root vegetables"           "canned beer"              
##  [5] "rolls/buns"                "whole milk"               
##  [7] "other vegetables"          "hygiene articles"         
##  [9] "whole milk"                "whole milk"               
## [11] "frozen fish"               "red/blush wine"           
## [13] "chocolate"                 "shopping bags"            
## [15] "dish cleaner"              "packaged fruit/vegetables"
## [17] "tropical fruit"            "rolls/buns"               
## [19] "cling film/bags"           "chocolate"                
## 
## $`1005`
## [1] "whipped/sour cream" "rolls/buns"         "margarine"         
## [4] "rolls/buns"

Con los datos agrupados ya en una lista, podemos convertirlos en una transacción, para esto, hacemos uso de la función as() del paquete arules:

datos.tra <- as(datos.lista, "transactions")
class(datos.tra)

## [1] "transactions"
## attr(,"package")
## [1] "arules"

summary(datos.tra)

## transactions as itemMatrix in sparse format with
##  3898 rows (elements/itemsets/transactions) and
##  167 columns (items) and a density of 0.05318557 
## 
## most frequent items:
##       whole milk other vegetables       rolls/buns             soda 
##             1783             1468             1360             1220 
##           yogurt          (Other) 
##             1100            27691 
## 
## element (itemset/transaction) length distribution:
## sizes
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   6 249  86 336 266 384 303 334 340 298 267 237 185 178 119  94  66  47  38  25 
##  21  22  23  24  25  26 
##  15  13   3   5   2   2 
## 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   6.000   8.000   8.882  12.000  26.000 
## 
## includes extended item information - examples:
##             labels
## 1 abrasive cleaner
## 2 artif. sweetener
## 3   baby cosmetics
## 
## includes extended transaction information - examples:
##   transactionID
## 1          1000
## 2          1001
## 3          1002

Según esto, tenemos 3898 transacciones y 167 ítems, la densidad es de 5.3% normal en estos casos. Ahora para visualizar la matriz de ítems usamos la función itemFrequencyPlot del paquete arules:

#Frecuencia Relativa
itemFrequencyPlot(datos.tra, topN = 5, type = "relative", horiz = TRUE, col = "purple")

#Frecuencia Absoluta
itemFrequencyPlot(datos.tra, topN = 5, type = "absolute", horiz = TRUE, col = "purple")

Para continuar viendo como quedo la lista de transacciones, podemos usar la función inspect del paquete arules, hay que aclarar que, si no se pone nada, va a mostrar todas las “facturas”, por lo que es mejor usar head, tail, o elegir una canasta aleatorea:

inspect(head(datos.tra, 2))

##     items                  transactionID
## [1] {canned beer,                       
##      hygiene articles,                  
##      misc. beverages,                   
##      pastry,                            
##      pickled vegetables,                
##      salty snack,                       
##      sausage,                           
##      semi-finished bread,               
##      soda,                              
##      whole milk,                        
##      yogurt}                        1000
## [2] {beef,                              
##      curd,                              
##      frankfurter,                       
##      rolls/buns,                        
##      sausage,                           
##      soda,                              
##      whipped/sour cream,                
##      white bread,                       
##      whole milk}                    1001

inspect(tail(datos.tra, 2))

##     items                    transactionID
## [1] {berries,                             
##      bottled water,                       
##      butter milk,                         
##      detergent,                           
##      herbs,                               
##      kitchen towels,                      
##      napkins,                             
##      newspapers,                          
##      onions,                              
##      other vegetables,                    
##      semi-finished bread,                 
##      tropical fruit,                      
##      whipped/sour cream,                  
##      yogurt}                          4999
## [2] {bottled beer,                        
##      fruit/vegetable juice,               
##      onions,                              
##      other vegetables,                    
##      root vegetables,                     
##      semi-finished bread,                 
##      soda}                            5000

inspect(datos.tra[102])

##     items               transactionID
## [1] {bottled water,                  
##      other vegetables,               
##      tropical fruit,                 
##      white bread,                    
##      whole milk,                     
##      yogurt}                     1105

Ahora, con el objeto creado, podemos empezar a crear las reglas de asociación, usando la función apriori del paquete arules y eligiendo los parámetros para hacer las reglas, estos parámetros son: umbral mínimo de soporte (supp), de confianza (conf), y el numero mínimo de la regla (minlen), que por defecto son 0.1, 0.8 y 1 respectivamente, pero para minlen es mejor elegir 2 porque tener 1 diría que aunque tenga la canasta vacía, aun así va a comprar el producto A. Antes de empezar, es ideal explicar bien que es cada cosa que nos dará la regla:

LHS: antecedente
RHS: consecuente
Support: frecuencia que aparece el itemset de la regla en el total de las transacciones
Confidence: probabilidad de comprar el RHS si se compro el LHS
Coverage: cuantas veces aparece el LHS como LHS (soporte del ítem LHS)
Lift: probabilidad que compren RHS con el LHS
Count: total de veces que se cumple la regla

reglas <- apriori(datos.tra, parameter = list(supp = 0.01, conf = 0.69, minlen = 2))

## Apriori
## 
## Parameter specification:
##  confidence minval smax arem  aval originalSupport maxtime support minlen
##        0.69    0.1    1 none FALSE            TRUE       5    0.01      2
##  maxlen target  ext
##      10  rules TRUE
## 
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
## 
## Absolute minimum support count: 38 
## 
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[167 item(s), 3898 transaction(s)] done [0.00s].
## sorting and recoding items ... [116 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 4 5 done [0.00s].
## writing ... [22 rule(s)] done [0.00s].
## creating S4 object  ... done [0.00s].

class(reglas)

## [1] "rules"
## attr(,"package")
## [1] "arules"

reglas

## set of 22 rules

Ahora, estas 22 reglas halladas en este caso no están organizadas de una manera específica, por lo que es necesario organizar las reglas antes de continuar, para esto organizamos las reglas por lift, por confianza y por soporte:

inspect(head(sort(reglas, by = "lift"), 5))

##     lhs                        rhs             support confidence   coverage     lift count
## [1] {domestic eggs,                                                                        
##      meat}                  => {whole milk} 0.01026167  0.7843137 0.01308363 1.714669    40
## [2] {chocolate,                                                                            
##      fruit/vegetable juice} => {whole milk} 0.01051821  0.7592593 0.01385326 1.659895    41
## [3] {bottled water,                                                                        
##      other vegetables,                                                                     
##      rolls/buns,                                                                           
##      yogurt}                => {whole milk} 0.01051821  0.7592593 0.01385326 1.659895    41
## [4] {bottled water,                                                                        
##      pip fruit,                                                                            
##      yogurt}                => {whole milk} 0.01026167  0.7407407 0.01385326 1.619410    40
## [5] {brown bread,                                                                          
##      rolls/buns,                                                                           
##      yogurt}                => {whole milk} 0.01282709  0.7352941 0.01744484 1.607502    50

inspect(head(sort(reglas, by = "confidence"), 5))

##     lhs                        rhs             support confidence   coverage     lift count
## [1] {domestic eggs,                                                                        
##      meat}                  => {whole milk} 0.01026167  0.7843137 0.01308363 1.714669    40
## [2] {chocolate,                                                                            
##      fruit/vegetable juice} => {whole milk} 0.01051821  0.7592593 0.01385326 1.659895    41
## [3] {bottled water,                                                                        
##      other vegetables,                                                                     
##      rolls/buns,                                                                           
##      yogurt}                => {whole milk} 0.01051821  0.7592593 0.01385326 1.659895    41
## [4] {bottled water,                                                                        
##      pip fruit,                                                                            
##      yogurt}                => {whole milk} 0.01026167  0.7407407 0.01385326 1.619410    40
## [5] {brown bread,                                                                          
##      rolls/buns,                                                                           
##      yogurt}                => {whole milk} 0.01282709  0.7352941 0.01744484 1.607502    50

inspect(head(sort(reglas, by = "support"), 5))

##     lhs                    rhs             support confidence   coverage     lift count
## [1] {butter,                                                                           
##      other vegetables,                                                                 
##      yogurt}            => {whole milk} 0.01564905  0.7093023 0.02206260 1.550679    61
## [2] {bottled beer,                                                                     
##      rolls/buns,                                                                       
##      yogurt}            => {whole milk} 0.01385326  0.7200000 0.01924064 1.574066    54
## [3] {other vegetables,                                                                 
##      rolls/buns,                                                                       
##      sausage,                                                                          
##      yogurt}            => {whole milk} 0.01359672  0.7066667 0.01924064 1.544917    53
## [4] {chocolate,                                                                        
##      pip fruit}         => {whole milk} 0.01334017  0.6933333 0.01924064 1.515767    52
## [5] {brown bread,                                                                      
##      rolls/buns,                                                                       
##      yogurt}            => {whole milk} 0.01282709  0.7352941 0.01744484 1.607502    50

Hemos encontrado 22 reglas, pero todavía no sabemos qué hacer con esto. Dentro de estas reglas existen reglas que son muy redundantes, es decir, que dicen que cosas que ya se pueden inferir por otra regla. Por esto, lo primero que hay que hacer es eliminar las reglas redundantes con las funciones is.redundant del paquete arules:

reglas.reduntantes <- is.redundant(reglas)
class(reglas.reduntantes)

## [1] "logical"

sum(reglas.reduntantes)

## [1] 0

reglas <- reglas[!reglas.reduntantes]

En este caso no habían reglas redundantes entonces todavía tenemos 22 reglas. Ahora, si se desea, se pueden convertir las reglas en un data frame usar dplyr:

reglas.df <- as(reglas, "data.frame")

Esto es útil si se desean graficar las reglas, como en este ejemplo:

plot(reglas, measure = c("support", "confidence"), shading = "lift")

Eliminamos el workspace.

Unidad 9: Clasificación

La tarea de clasificación tiene como finalidad de crear categorías por medio de una muestra a individuos u objetos con la intención de que cuando entre otro individuo a la muestra este caiga en una de las categorías ya creadas. En esta unidad nos enfocaremos en solo dos modelos:

Modelo Logit
Modelo kNN (k-Near Neighbors)

Para ambos modelos es necesario destinar una parte de la muestra para la estimación del modelo y otra para la evaluación (generalmente se divide 80/20) para evitar el overfitting del modelo. Otra cosa importante para el modelo es el indice de validación cruzada pues este nos dirá que tan bueno es el modelo antes del evaluarlo.

Antes de empezar utilizaremos dos bases para este tema por dos razones, la de logit tiene variables de dummy que ayudara a explicar eso para que sirve y la otra la usaremos para el kNN porque el kNN se tarda en cargar, de igual manera vamos a crear los indices para la estimación, de evaluación y de validación cruzada:

require(pacman)
p_load(tidyverse, fastDummies, ROCit, descr, caret, mfx, rio)

base_logit <- import("bases/bank-additional-full.csv")
base_kNN <- import("bases/data_train.RData")

#Estimación y evaluación del logit
set.seed(171124)

est_index_logit <- sample(1:nrow(base_logit), 0.8 * nrow(base_logit))
eval_index_logit <- setdiff(1:nrow(base_logit), est_index_logit)

#Estimación y evaluación ara el kNN
set.seed(171124)

est_index_kNN <- sample(1:nrow(base_kNN), 0.8 * nrow(base_kNN))
eval_index_kNN <- setdiff(1:nrow(base_kNN), est_index_kNN)

#Validación cruzada
set.seed(171124)

cross_validation <- trainControl(method = "cv", number = 5)

Ya con todo esto definido. Podemos empezar con el modelo logit. Este modelo consiste en una generalización del modelo de regresión lineal múltiple. Esto implica que hay una variable cuantitativa continua, pero en este caso, la variable sera cualitativa.

Antes de empezar bien con el análisis, eliminaremos 3 variables que no nos aportaran con el análisis del modelo:

base_logit <- subset(base_logit, select = c(-duration, -pdays, -nr.employed))
glimpse(base_logit)

## Rows: 41,188
## Columns: 18
## $ age            <int> 56, 57, 37, 40, 56, 45, 59, 41, 24, 25, 41, 25, 29, 57,…
## $ job            <chr> "housemaid", "services", "services", "admin.", "service…
## $ marital        <chr> "married", "married", "married", "married", "married", …
## $ education      <chr> "basic.4y", "high.school", "high.school", "basic.6y", "…
## $ default        <chr> "no", "unknown", "no", "no", "no", "unknown", "no", "un…
## $ housing        <chr> "no", "no", "yes", "no", "no", "no", "no", "no", "yes",…
## $ loan           <chr> "no", "no", "no", "no", "yes", "no", "no", "no", "no", …
## $ contact        <chr> "telephone", "telephone", "telephone", "telephone", "te…
## $ month          <chr> "may", "may", "may", "may", "may", "may", "may", "may",…
## $ day_of_week    <chr> "mon", "mon", "mon", "mon", "mon", "mon", "mon", "mon",…
## $ campaign       <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ previous       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ poutcome       <chr> "nonexistent", "nonexistent", "nonexistent", "nonexiste…
## $ emp.var.rate   <dbl> 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, …
## $ cons.price.idx <dbl> 93.994, 93.994, 93.994, 93.994, 93.994, 93.994, 93.994,…
## $ cons.conf.idx  <dbl> -36.4, -36.4, -36.4, -36.4, -36.4, -36.4, -36.4, -36.4,…
## $ euribor3m      <dbl> 4.857, 4.857, 4.857, 4.857, 4.857, 4.857, 4.857, 4.857,…
## $ y              <chr> "no", "no", "no", "no", "no", "no", "no", "no", "no", "…

Ahora con la base ya limpiada, es necesario convertir la variables que no son cuantitativas en dummies para que su respuesta sean solamente 0 y 1 y se puedan así interpretar como unas variables numéricas, el proceso entero es el siguiente:

dummies_job <- dummy_cols(base_logit$job, remove_first_dummy = TRUE, 
                          remove_selected_columns = TRUE, omit_colname_prefix = TRUE)
names(dummies_job) <- c("j_blue_collar", "j_entrepreneur", "j_housemaid", "j_management",
                        "j_retired", "j_self_employed", "j_services", "j_student", 
                        "j_technician", "j_unemployed", "j_unknown")

dummies_marital <- dummy_cols(base_logit$marital, remove_first_dummy = TRUE, 
                              remove_selected_columns = TRUE)
names(dummies_marital) <- c("m_married", "m_single", "m_unknown")

dummies_education <- dummy_cols(base_logit$education, remove_first_dummy = TRUE, 
                                remove_selected_columns = TRUE)
names(dummies_education) <- c("e_basic_6y", "e_basic_9y", "e_high_school", "e_illiterate",
                              "e_professional_course", "e_university_degree", "e_unknown")

dummies_default <- dummy_cols(base_logit$default, remove_first_dummy = TRUE,
                              remove_selected_columns = TRUE)
names(dummies_default) <- c("def_unknown", "def_yes")

dummies_housing <- dummy_cols(base_logit$housing, remove_first_dummy = TRUE,
                              remove_selected_columns = TRUE)
names(dummies_housing) <- c("h_unknown", "h_yes")

dummies_loan <- dummy_cols(base_logit$loan, remove_first_dummy = TRUE,
                           remove_selected_columns = TRUE)
names(dummies_loan) <- c("l_unknown", "l_yes")

dummies_contact <- dummy_cols(base_logit$contact, remove_first_dummy = TRUE,
                              remove_selected_columns = TRUE)
names(dummies_contact) <- c("c_telephone")

dummies_month <- dummy_cols(base_logit$month, remove_first_dummy = TRUE,
                            remove_selected_columns = TRUE)
names(dummies_month) <- c("mon_aug", "mon_dec", "mon_jul", "mon_jun",
                          "mon_mar", "mon_may", "mon_nov", "mon_oct", "mon_sep")

dummies_dayWeek <- dummy_cols(base_logit$day_of_week, remove_first_dummy = TRUE,
                              remove_selected_columns = TRUE)
names(dummies_dayWeek) <- c("d_mon", "d_thu", "d_tue", "d_wed")

datos_dummies <- subset(base_logit, select = -c(job, marital, education, default, housing,
                                                loan, contact, month, day_of_week))

datos_dummies <- cbind(datos_dummies, dummies_job, dummies_marital, dummies_education, 
                       dummies_default, dummies_housing, dummies_loan, dummies_contact, 
                       dummies_month, dummies_dayWeek)

datos_dummies$y <- relevel(as.factor(datos_dummies$y), ref = "no")

Como ya tenemos las bases con las variables dummy podemos pasar a sacar las bases de estimación y evaluación:

set.seed(171124)

datos_dummies_est <- datos_dummies[est_index_logit, ]
datos_dummies_eval <- datos_dummies[eval_index_logit, ]

Con todo esto ya listo, podemos empezar a realizar el modelo, para este tenemos el modelo máximo que tiene todas la variables y el mínimo que no tiene ninguna, lo que hará el algoritmo es añadir variables hasta que llegue al punto optimo:

max_modelo <- glm(y ~ ., data = datos_dummies_est, family = "binomial")
min_modelo <- glm(y ~ 1, data = datos_dummies_est, family = "binomial")

modelo_forward <- step(min_modelo, 
                              scope = list(lower = min_modelo,
                                           upper = max_modelo),
                              direction = "forward")

modelo_logit <- glm(modelo_forward$formula, data = datos_dummies_est, family = "binomial")

Pero que hacemos con esto? Bueno, hay varios indices que ayudan evaluar el modelo. El primero es el LRI, que es un indice de razón de verosimilitud, entre mayor sea el indice (entre 0 y 1) esta más ajustado el modelo.

1 - logLik(modelo_logit)/logLik(min_modelo)

## 'log Lik.' 0.2158803 (df=26)

El segundo indice que se puede utilizar es la R cuadrada que es muy similar al LRI

LogRegR2(modelo_logit)

## Chi2                 5006.4 
## Df                   25 
## Sig.                 0 
## Cox and Snell Index  0.1409596 
## Nagelkerke Index     0.2789597 
## McFadden's R2        0.2158803

Ahora, una diferencia del modelo logit y del kNN es que este no te da los valores directamente, sino la probabilidad de que sea ese valor por lo que hay que sacar las probabilidades y por medio de la curva ROC y el indice AUC se pueden conseguir esas predicciones reales:

prob_pred_modelo <- predict(modelo_logit, newdata = datos_dummies_eval, type = "response")

La curva ROC es una herramienta para evaluar el desempeño de un modelo de clasificación, el eje x corresponde a la tasa de falsos positivos, el eje y es la tasa de verdaderos positivos y la resta es un modelo aleatorio, es decir, entre más lejos de la recta y más cerca de la esquina superior izquierda mejor:

roc_modelo <- rocit(score = prob_pred_modelo, class = datos_dummies_eval$y, method = "bin")
plot(roc_modelo, values = TRUE)

summary(roc_modelo)

##                             
##  Method used: binormal      
##  Number of positive(s): 930 
##  Number of negative(s): 7308
##  Area under curve: 0.7843

El área bajo de la curva es el AUC, este me dice que tan bien predice si el cliente dice “si” o dice “no”. Ahora, vamos a ver cual es valor de corte optimo, que se puede visualizar en el plot, pero no sabemos ese valor, este valor nos dirá la probabilidad mínima de que el cliente diga que “si”.

metricas_modelo <- measureit(score = modelo_logit$fitted.values, class = datos_dummies_est$y,
                              measure = c("SENS", "SPEC", "PREC", "NPV", "ACC", "FSCR", "FPR"))

opt_pos_modelo <- which.max(metricas_modelo$SENS + metricas_modelo$SPEC - 1)
v_corte_modelo <- metricas_modelo$Cutoff[opt_pos_modelo]
v_corte_modelo

## [1] 0.1184166

Con esto ya podemos construir la métricas con el mejor punto de corte y construir la matriz de confusión que nos dará el numero de casos de que el dato sea un falso positivo, un falso negativo, un verdadero positivo o un verdadero negativo:

pred_modelo <- ifelse(prob_pred_modelo > v_corte_modelo, "yes", "no") 
pred_modelo <- as.factor(pred_modelo)

mo_logit <- confusionMatrix(pred_modelo, datos_dummies_eval$y, positive = "yes")

metrica_modelo <- cbind(mo_logit$byClass[1], 
                         mo_logit$byClass[3], 
                         mo_logit$byClass[5], 
                         mo_logit$byClass[4], 
                         mo_logit$overall[1])

tabla <- metrica_modelo
tabla <- as.data.frame(tabla)
names(tabla) <- c("SENS", "SPEC", "PREC", "NPV", "ACC")

mo_logit$table

##           Reference
## Prediction   no  yes
##        no  6211  346
##        yes 1097  584

glimpse(tabla)

## Rows: 1
## Columns: 5
## $ SENS <dbl> 0.627957
## $ SPEC <dbl> 0.3474123
## $ PREC <dbl> 0.3474123
## $ NPV  <dbl> 0.947232
## $ ACC  <dbl> 0.8248361

Sensibilidad: El modelo tiene una excelente capacidad para identificar correctamente los casos positivos. Detecta la gran mayoría de ellos.
Especificidad: Su habilidad para reconocer los casos negativos es limitada, aunque mejora respecto a otros modelos con especificidad muy baja.
Precisión: De todos los casos que el modelo clasifica como positivos, menos de la mitad lo son realmente. Esto indica una proporción considerable de falsos positivos.
Valor predictivo negativo: El modelo es altamente confiable cuando predice que un caso es negativo. Casi todos los negativos predichos realmente lo son.
Exactitud: El porcentaje total de clasificaciones correctas es bueno. El modelo acierta en una proporción significativa de las predicciones.

Ahora bien, podemos ver la relación entre las variables utilizando sus valores marginales, para esto podemos usar la librería mfx.

efectos_marginales <- logitmfx(modelo_logit$formula, 
        data = datos_dummies_est, atmean = F)
print(efectos_marginales)

## Call:
## logitmfx(formula = modelo_logit$formula, data = datos_dummies_est, 
##     atmean = F)
## 
## Marginal Effects:
##                           dF/dx   Std. Err.        z     P>|z|    
## euribor3m            0.02093722  0.00635755   3.2933 0.0009902 ***
## mon_may             -0.03475553  0.00479847  -7.2430 4.387e-13 ***
## poutcomenonexistent  0.03372135  0.00605125   5.5726 2.509e-08 ***
## poutcomesuccess      0.20071480  0.01517277  13.2286 < 2.2e-16 ***
## cons.price.idx       0.11178116  0.00834113  13.4012 < 2.2e-16 ***
## cons.conf.idx        0.00163694  0.00038708   4.2289 2.348e-05 ***
## mon_mar              0.13860805  0.01667646   8.3116 < 2.2e-16 ***
## c_telephone         -0.04413479  0.00484371  -9.1118 < 2.2e-16 ***
## emp.var.rate        -0.08805695  0.00857463 -10.2695 < 2.2e-16 ***
## d_mon               -0.01813065  0.00374508  -4.8412 1.291e-06 ***
## mon_nov             -0.02984090  0.00582421  -5.1236 2.998e-07 ***
## mon_jun             -0.02456177  0.00611106  -4.0192 5.839e-05 ***
## def_unknown         -0.01756662  0.00452448  -3.8826 0.0001034 ***
## campaign            -0.00335804  0.00079839  -4.2060 2.599e-05 ***
## j_retired            0.02926092  0.00754104   3.8802 0.0001044 ***
## mon_aug              0.02992973  0.00837759   3.5726 0.0003535 ***
## j_student            0.01882848  0.00917330   2.0525 0.0401181 *  
## d_wed                0.00973390  0.00411475   2.3656 0.0180004 *  
## mon_jul              0.01610475  0.00673655   2.3907 0.0168185 *  
## j_blue_collar       -0.00679688  0.00440031  -1.5446 0.1224344    
## h_unknown           -0.01906381  0.00960850  -1.9841 0.0472496 *  
## mon_dec              0.02988029  0.01805847   1.6546 0.0979972 .  
## m_single             0.00559765  0.00358680   1.5606 0.1186119    
## previous             0.00687521  0.00467744   1.4699 0.1415982    
## e_university_degree  0.00524954  0.00361056   1.4539 0.1459627    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
##  [1] "mon_may"             "poutcomenonexistent" "poutcomesuccess"    
##  [4] "mon_mar"             "c_telephone"         "d_mon"              
##  [7] "mon_nov"             "mon_jun"             "def_unknown"        
## [10] "j_retired"           "mon_aug"             "j_student"          
## [13] "d_wed"               "mon_jul"             "j_blue_collar"      
## [16] "h_unknown"           "mon_dec"             "m_single"           
## [19] "e_university_degree"

Finalmente podemos realizar una validación cruzada del modelo logit para ver que tan bueno es. Esto se hace incluso antes de la matriz de confusión y se realiza de la siguiente manera:

modelo_logit_cv <- train(form = modelo_logit$formula, data = datos_dummies_est, 
                         method = "glm", family = "binomial", trControl = cross_validation)
modelo_logit_cv

## Generalized Linear Model 
## 
## 32950 samples
##    24 predictor
##     2 classes: 'no', 'yes' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 26360, 26360, 26360, 26360, 26360 
## Resampling results:
## 
##   Accuracy  Kappa    
##   0.899393  0.2964151

Esto nos dara el accuracy del modelo evaluado consigo mismo y el kappa (concordancia entre la predicción y los valores reales) de la misma manera.

Ahora para el kNN tenemos que es modelo es un algoritmo de aprendizaje supervisado que se basa en la proximidad de los datos. Entonces lo primero que debemos hacer es estandarizar las variables (claramente solo las cuantitativas):

datos_expl <- base_kNN %>% 
  dplyr::select(-label, -id_row) %>%
  scale()
head(datos_expl, 3)

##    num_words num_unique_words num_stopwords  num_links num_unique_domains
## 1 -0.2209053       -0.2844703    -0.1711331  0.2619614          1.1495860
## 2 -0.2730243       -0.3731394    -0.3232038 -0.1503332         -0.2269033
## 3  1.1007776        2.3534350     1.0285363 -0.1503332         -0.2269033
##   num_email_addresses num_spelling_errors num_urgent_keywords
## 1          -0.1530115          -0.3038335          -0.4565362
## 2          -0.1530115          -0.3699350          -0.4565362
## 3          -0.1530115           0.5720118           2.8045283

base_kNN$label <- factor(base_kNN$label, levels = c(0, 1), labels = c("Seguro", "Phishing"))

Como hicimos en el modelo Logit, vamos a destinar el 80% de la muestra para entrenar el algoritmo (la estimación) y el otro 20% para la evaluación:

datos_expl_est <- datos_expl[est_index_kNN,]
label_est <- base_kNN$label[est_index_kNN]

datos_expl_eval <- datos_expl[eval_index_kNN,]
label_eval <- base_kNN$label[eval_index_kNN]
datos_expl_eval <- as.data.frame(datos_expl_eval)

Ahora podemos correr el algoritmo KNN, y para mejorar la precisión de este lo corremos realizando una validación cruzada:

set.seed(171124)

kNN_res <- train(datos_expl_est, label_est, method = "knn", preProcess = NULL,
                 tuneGrid = expand.grid(k = seq(3, 10, 1)), trControl = cross_validation)
kNN_res$bestTune$k

## [1] 7

Con esto ya se pueden hacer la predicciones y la matriz de confusión, porque no hay que hacer lo que se hace en el modelo logit (ya se menciono esto previamente):

pred_kNN <- predict(kNN_res, newdata = datos_expl_eval, type = "raw")

mc_kNN <- confusionMatrix(pred_kNN, label_eval, positive = "Phishing")
mc_kNN$table

##           Reference
## Prediction Seguro Phishing
##   Seguro     1413      364
##   Phishing    385      889

kNN <- cbind(mc_kNN$byClass[1], mc_kNN$byClass[3], mc_kNN$byClass[5], 
               mc_kNN$byClass[4], mc_kNN$overall[1])

mc_kNN$overall

##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##   7.545067e-01   4.941153e-01   7.388318e-01   7.696922e-01   5.893150e-01 
## AccuracyPValue  McnemarPValue 
##   6.840787e-82   4.649110e-01

mc_kNN$byClass

##          Sensitivity          Specificity       Pos Pred Value 
##            0.7094972            0.7858732            0.6978022 
##       Neg Pred Value            Precision               Recall 
##            0.7951604            0.6978022            0.7094972 
##                   F1           Prevalence       Detection Rate 
##            0.7036011            0.4106850            0.2913799 
## Detection Prevalence    Balanced Accuracy 
##            0.4175680            0.7476852

kNN <- cbind(mc_kNN$overall[1], mc_kNN$byClass[1], mc_kNN$byClass[2], 
               mc_kNN$byClass[5], mc_kNN$byClass[7], mc_kNN$overall[2])

kNN <- as.data.frame(kNN)
names(kNN) <- c("ACC", "SENS", "SPEC", "PREC", "F1", "KAPPA")
row.names(kNN) <- NULL
kNN

##         ACC      SENS      SPEC      PREC        F1     KAPPA
## 1 0.7545067 0.7094972 0.7858732 0.6978022 0.7036011 0.4941153

Borramos workspace.

Unidad 10: Detección de Anomalías

La ultima unidad del curso se enfoca en la detección de anomalías. Esta tarea se enfoca es como su nombre lo implica, por medio de un análisis exploratorio de los datos, detectar anomalías como lo pueden ser los outliers. En este curso nos enfocamos en las anomalías univariadas, estas son anomalías que ocurren en una sola variable. Para hallar estas anomalías utilizaremos un EDA (exploratory data analysis) para:

Visualizar los datos, para entender gráficamente el comportamiento de la muestra
calcular estadísticas, como estadísticas de tendencia central, dispersión, posición y forma

De igual manera veremos sobre la ley de Benford para la detección de fraudes.

Antes de empezar el EDA y la detección de fraudes, necesitamos cargar los paquetes y librerias necesarios:

require(pacman)
p_load(tidyverse, robustbase, grDevices, funModeling, univOutl, qqplotr,
       gridExtra, outliers, EnvStats, benford.analysis, BenfordTests, rio)

datos23 <- import("bases/datos_credito.RData")
glimpse(datos23)

## Rows: 1,000
## Columns: 14
## $ cuenta.corriente      <fct> less than 0 DM, 0 or less than 200 DM, no checki…
## $ historial.crediticio  <fct> critical account/ other credits existing (not at…
## $ proposito             <fct> radio/television, radio/television, education, f…
## $ monto.credito         <int> 1169, 5951, 2096, 7882, 4870, 9055, 2835, 6948, …
## $ cuenta.ahorros        <fct> unknown/ no savings account, less than 100 DM, l…
## $ porcentaje.disponible <int> 4, 2, 2, 2, 3, 2, 3, 2, 2, 4, 3, 3, 1, 4, 2, 4, …
## $ estadocivil.sexo      <fct> male : single, female : divorced/separated/marri…
## $ residente             <int> 4, 2, 3, 4, 4, 4, 4, 2, 4, 2, 1, 4, 1, 4, 4, 2, …
## $ edad                  <int> 67, 22, 49, 45, 53, 35, 53, 35, 61, 28, 25, 24, …
## $ propiedad.hogar       <fct> own, own, own, for free, for free, for free, own…
## $ numero.creditos       <int> 2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, …
## $ ocupacion             <fct> skilled employee / official, skilled employee / …
## $ dependientes          <int> 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ trabajador.extranjero <fct> yes, yes, yes, yes, yes, yes, yes, yes, yes, yes…

datos6 <- import("bases/datosBenford.txt")
glimpse(datos6)

## Rows: 184,412
## Columns: 4
## $ registro  <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ VendorNum <int> 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, …
## $ Date      <chr> "2/01/10", "2/01/10", "2/01/10", "2/01/10", "2/01/10", "2/01…
## $ Amount    <dbl> 36.08, 77.80, 34.97, 59.00, 59.56, 50.38, 26.57, 102.17, 25.…

Como ya dijimos, el primer paso es la visualización de los datos. Para esto utilizaremos la variable de monto_credito de la base para los capítulos 2 y 3. Las visualizaciones principales para este proceso son:

Histograma
Boxplot
Boxplot ajustado (toma en cuenta la asimetría de los datos)

#Histograma
ggplot(datos23, aes(x = monto.credito)) +
  geom_histogram(bins = round(sqrt(length(datos23$monto.credito))), 
                 fill = "purple", col = "black", alpha = 0.7) +
  labs(title = "Distribución de los Montos de los Creditos", x = "Monto de Credito", 
       y = "Numero de Clientes", caption = "Fuente: Elaboración Propia") +
  theme_bw()

#Boxplot
ggplot(datos23, aes(x = monto.credito)) +
  geom_boxplot(fill = "purple", col = "black", alpha = 0.5, outlier.color = "black") +
  labs(title = "Distribución de los Montos de los Creditos", x = "Monto de Credito", 
       y = "Numero de Clientes", caption = "Fuente: Elaboración Propia") +
  theme_bw()

#Boxplot Ajustado
adjbox_stats <- adjboxStats(datos23$monto.credito)$stats
ggplot(datos23, aes(x = monto.credito, y = "")) +
  geom_boxplot(xmin = adjbox_stats[1], 
               xmax = adjbox_stats[5], middle = adjbox_stats[3],
               upper = adjbox_stats[4], lower = adjbox_stats[2], outlier.shape = NA,
               fill = "purple", col = "black", alpha = 0.5) +
  geom_point(data = subset(datos23, monto.credito < adjbox_stats[1] | 
                             monto.credito > adjbox_stats[5]),
             size = 2, alpha = 0.5) +
  labs(title = "Distribución Ajustada de los Montos de los Creditos", x = "Monto de Credito", 
       y = "Numero de Clientes", caption = "Fuente: Elaboración Propia") +
  theme_bw()

Esto es todo por el análisis gráfico. Esto nos muestra que incluso tomando el cuenta el sesgo de los datos, existen datos atípicos, por lo que pasamos a analizar las métricas para detectar outliers. Las métricas para detectar esto son:

Uso de Percentiles: el uso de los percentiles como umbrales para la detección de outliers, generalmente se utilizan los percentiles 2.5 y 97.5 pero se pueden usar otros como 99 y 1. El problema con este que siempre va a hallar atípicos
Z-Score: es un indicador de la desviación estándar por lo que se enfoca en estandarizar los datos y según esto entre más lejano de 0 este el valor, más inusual es. Normalmente se utilizan de umbrales a -3 y a 3, sin embargo, si los datos no se distribuyen de manera normal no puede ser conveniente
Rango Intercuartilico: es la prueba de tukey utilizada en la unidad 4. Es similar al boxplot tradicional y ajustado que se utilizaron en el análisis gráfico
Método de Hampel: por medio del uso de la desviación absoluta de la mediana (MAD) se hallan outliers según la distancia que tiene con sus vecinos, aquellos puntos que se desvíen más de lo normal, son atípicos. Todo depende de la K que multiplica la MAD, si es muy pequeña encontrara muchos y si es muy grande no encontrara.
Métricas que Emplean Estimadores Robustos de Varianza: estos estimadores hacen que los atípicos no impacten tanto en el resultado de la media y de la varianza, lo que permite hacer un proceso de escalar los datos, similar al Z-score, para poderlos comparar a una distribución estadística

#Percentil
#datos encima del 99.5% del resto
datos23 %>%
  dplyr::select(monto.credito) %>%
  filter(monto.credito > quantile(monto.credito, 0.995))

##   monto.credito
## 1         15945
## 2         15653
## 3         15857
## 4         15672
## 5         18424

#Z-Score
z_scores <- as.data.frame(scale(datos23$monto.credito))
names(z_scores) <- "score"

z_scores %>%
  filter(score < -3) %>%
  dim() #datos atipicamente bajos, el dim dice solo la cantidad de observaciones

## [1] 0 1

z_scores %>%
  filter(score > 3) #datos atipicamente altos, sin dim da los z_scores

##       score
## 1  3.297418
## 2  3.949976
## 3  3.309108
## 4  4.489877
## 5  3.070333
## 6  3.997447
## 7  3.152168
## 8  3.091589
## 9  3.714389
## 10 4.077866
## 11 3.913486
## 12 3.438061
## 13 3.007274
## 14 3.230107
## 15 3.164568
## 16 4.386432
## 17 3.810395
## 18 3.864243
## 19 3.333198
## 20 4.458702
## 21 3.027112
## 22 4.393163
## 23 5.368103
## 24 4.118252
## 25 3.357643

#IQR
boxB(datos23$monto.credito, method = "resistant") #manera tradicional (out)

## No. of outliers in left tail: 0

## No. of outliers in right tail: 72

## $quartiles
##     25%     50%     75% 
## 1365.50 2319.50 3972.25 
## 
## $fences
##     lower     upper 
## -2544.625  7882.375 
## 
## $excluded
## integer(0)
## 
## $outliers
##  [1]   6  18  19  58  64  71  79  88  96 106 131 135 137 181 206 227 237 269 273
## [20] 275 286 292 296 305 334 374 375 379 382 396 403 418 432 451 492 497 510 526
## [39] 550 564 616 617 638 646 654 658 673 685 715 737 745 764 772 806 809 813 819
## [58] 829 833 855 882 888 896 903 916 918 922 928 946 954 981 984
## 
## $lowOutl
## integer(0)
## 
## $upOutl
##  [1]   6  18  19  58  64  71  79  88  96 106 131 135 137 181 206 227 237 269 273
## [20] 275 286 292 296 305 334 374 375 379 382 396 403 418 432 451 492 497 510 526
## [39] 550 564 616 617 638 646 654 658 673 685 715 737 745 764 772 806 809 813 819
## [58] 829 833 855 882 888 896 903 916 918 922 928 946 954 981 984

boxB(datos23$monto.credito, method = "adjbox") #manera ajustada (out_adj)

## Warning in boxB(datos23$monto.credito, method = "adjbox"): With method='adjbox'
## the argument k is set equal to 1.5

## The MedCouple skewness measure is: 0.3771

## No. of outliers in left tail: 18

## No. of outliers in right tail: 1

## $quartiles
##     25%     50%     75% 
## 1365.50 2319.50 3972.25 
## 
## $fences
##      lower      upper 
##   499.6755 16097.1160 
## 
## $excluded
## integer(0)
## 
## $outliers
##  [1]  27  28  40 112 158 178 250 310 380 459 472 494 591 722 726 751 812 916 965
## 
## $lowOutl
##  [1]  27  28  40 112 158 178 250 310 380 459 472 494 591 722 726 751 812 965
## 
## $upOutl
## [1] 916

LocScaleB(datos23$monto.credito, method = "IQR") #otra manera de ajuste

## No. of outliers in left tail: 0

## No. of outliers in right tail: 67

## $pars
##   median    scale 
## 2319.500 1932.357 
## 
## $bounds
## lower.low  upper.up 
## -3477.572  8116.572 
## 
## $excluded
## integer(0)
## 
## $outliers
##  [1]   6  19  58  64  71  79  88  96 106 131 135 137 181 206 227 237 269 273 275
## [20] 286 292 296 305 334 374 375 379 382 396 403 418 432 451 492 497 510 550 564
## [39] 616 617 638 658 673 685 715 737 745 764 806 809 813 819 829 833 855 882 888
## [58] 896 903 916 918 922 928 946 954 981 984
## 
## $lowOutl
## integer(0)
## 
## $upOutl
##  [1]   6  19  58  64  71  79  88  96 106 131 135 137 181 206 227 237 269 273 275
## [20] 286 292 296 305 334 374 375 379 382 396 403 418 432 451 492 497 510 550 564
## [39] 616 617 638 658 673 685 715 737 745 764 806 809 813 819 829 833 855 882 888
## [58] 896 903 916 918 922 928 946 954 981 984

#MAD
LocScaleB(datos23$monto.credito, method = "MAD")

## No. of outliers in left tail: 0

## No. of outliers in right tail: 100

## $pars
##   median    scale 
## 2319.500 1627.153 
## 
## $bounds
## lower.low  upper.up 
##  -2561.96   7200.96 
## 
## $excluded
## integer(0)
## 
## $outliers
##   [1]   4   6  18  19  49  58  64  71  79  88  96 106 109 114 131 135 137 154
##  [19] 164 176 181 206 227 228 237 256 269 273 275 286 288 292 295 296 305 333
##  [37] 334 374 375 376 379 382 388 396 403 412 418 432 451 468 492 497 510 526
##  [55] 539 550 564 616 617 638 646 651 654 658 673 685 715 716 737 745 764 772
##  [73] 797 805 806 809 813 816 819 829 833 855 869 871 881 882 888 890 896 903
##  [91] 916 918 922 928 946 954 972 974 981 984
## 
## $lowOutl
## integer(0)
## 
## $upOutl
##   [1]   4   6  18  19  49  58  64  71  79  88  96 106 109 114 131 135 137 154
##  [19] 164 176 181 206 227 228 237 256 269 273 275 286 288 292 295 296 305 333
##  [37] 334 374 375 376 379 382 388 396 403 412 418 432 451 468 492 497 510 526
##  [55] 539 550 564 616 617 638 646 651 654 658 673 685 715 716 737 745 764 772
##  [73] 797 805 806 809 813 816 819 829 833 855 869 871 881 882 888 890 896 903
##  [91] 916 918 922 928 946 954 972 974 981 984

#RRMS
out_RRMS <- hotspots::outliers(datos23$monto.credito, tail = "both")
out_RRMS #umbrales

## outlier cutoff (positive): 
## [1] 6711.865
## outlier cutoff (negative): 
## [1] -2072.865

summary(out_RRMS) #podemos ver los outliers superiores e inferiores

## 
## Source data: datos23$monto.credito
## Distribution and probability: t, 0.99
## Tail: positive and negative outliers                           
## Mean:             3.271e+03
## Median:           2.320e+03
## Min:              2.500e+02
## Max:              1.842e+04
## mad:              1.627e+03
## CV (mad/median):  7.015e-01
## 
## n =  1000
## 
## positive outliers:
##   Cutoff  number positive outliers  % positive outliers   % sum
##     6712                       117                  11.7   34.1
## 
## 
## negative outliers:
##   Cutoff  number negative outliers  % negativeoutliers   % sum
##    -2073                         0                    0      0

plot(out_RRMS)

out_ind_RRMS <- which(datos23$monto.credito < out_RRMS$negative.cut 
                      | datos23$monto.credito > out_RRMS$positive.cut)
out_ind_RRMS

##   [1]   4   6   8  18  19  30  49  58  64  71  79  88  96 100 106 109 114 117
##  [19] 131 132 135 137 154 155 164 176 181 206 227 228 237 256 269 273 275 286
##  [37] 288 292 295 296 305 333 334 374 375 376 379 382 388 396 403 412 418 432
##  [55] 451 468 492 497 508 510 518 523 526 539 550 553 564 570 616 617 638 646
##  [73] 651 654 658 673 685 715 716 737 739 745 764 772 797 805 806 809 813 816
##  [91] 819 829 833 847 855 869 871 880 881 882 888 890 896 903 916 918 922 925
## [109] 928 940 946 954 969 972 974 981 984

out_RRMS_dataset <- datos23[out_ind_RRMS, ]

Ya hemos hablado de como reconocer posibles outliers, pero existen unas pruebas estadísticas que también se encargan de detectar anomalías, sin embargo es necesario que estas variables cumplan con el supuesto de normalidad, para verificar esto, podemos hacer una visualización donde comparamos una distribución con los datos de la variable, si los datos están en la linea o al menos dentro del intervalo de confianza, entonces se puede realizar, si no, no.

ggplot(datos23, aes(sample = edad)) +
  stat_qq_band() +
  stat_qq_point(color = "purple", alpha = 0.5) +
  stat_qq_line() +
  labs(title = "Comprobación de Supuesto de Normalidad para Edad", 
       x = "Cuantiles Teóricos",
       y = "Cuantiles Observados", caption = "Fuente: Elaboración Propia") +
  theme_bw() #cerca a una distribución normal

ggplot(datos23, aes(sample = porcentaje.disponible)) +
  stat_qq_band() +
  stat_qq_point(color = "purple", alpha = 0.5) +
  stat_qq_line() +
  labs(title = "Comprobación de Supuesto de Normalidad para Porcentaje Disponible", 
       x = "Cuantiles Teóricos", y = "Cuantiles Observados", 
       caption = "Fuente: Elaboración Propia") +
  theme_bw() #no es normal

ggplot(datos23, aes(sample = residente)) +
  stat_qq_band() +
  stat_qq_point(color = "purple", alpha = 0.5) +
  stat_qq_line() +
  labs(title = "Comprobación de Supuesto de Normalidad para Residente", 
       x = "Cuantiles Teóricos", y = "Cuantiles Observados", 
       caption = "Fuente: Elaboración Propia") +
  theme_bw() #no es normal

ggplot(datos23, aes(sample = numero.creditos)) +
  stat_qq_band() +
  stat_qq_point(color = "purple", alpha = 0.5) +
  stat_qq_line() +
  labs(title = "Comprobación de Supuesto de Normalidad para Numero de Creditos", 
       x = "Cuantiles Teóricos", y = "Cuantiles Observados", 
       caption = "Fuente: Elaboración Propia") +
  theme_bw() #no es normal

Como ninguno tiene una distribución normal, entonces crearemos nosotros los datos, con unos atipicos intencionales y realizaremos las pruebas:

set.seed(171124)

datos_prueba <- data.frame(W = c(rnorm(999), 6), Z = c(rnorm(993), -9, -6, -10, 10, 8, 6, 7.5))

p1 <- ggplot(datos_prueba, aes(sample = W)) +
  stat_qq_band() +
  stat_qq_point(color = "purple", alpha = 0.5, size = 2) +
  stat_qq_line() +
  labs(title = "Comprobación de Supuesto de Normalidad", x = "Cuantiles Teóricos",
       y = "Cuantiles Observados", caption = "Fuente: Elaboración Propia") +
  theme_bw()

p2 <- ggplot(datos_prueba, aes(sample = Z)) +
  stat_qq_band() +
  stat_qq_point(color = "purple", alpha = 0.5, size = 2) +
  stat_qq_line() +
  labs(title = "Comprobación de Supuesto de Normalidad", x = "Cuantiles Teóricos",
       y = "Cuantiles Observados", caption = "Fuente: Elaboración Propia") +
  theme_bw()

grid.arrange(p1, p2, ncol = 2)

Con los datos ya obtenidos podemos realizar las siguientes pruebas:

Prueba de Grubbs: es una prueba empleada para detectar un único outlier. Esta prueba sigue las siguientes hipótesis, la nula dice que la muestra proviene de una distribución normal sin outliers y la alterna dice que hay un outlier en la muestra. Este test funciona con muestras grandes así no sean normales, pero si son pequeñas va a ser muy sensible. Si se intuye la existencia de más de un outlier, esta prueba no sirve.
Prueba de Dixon: hace lo mismo que Grubbs solo que recibe muestras entre 3 y 25
Prueba de Rosner: aunque tiene el mismo propósito de que la prueba de Grubbs y Dixon, esta si puede detectar varios outliers.

#Grubbs
grubbs.test(datos_prueba$W, type = 10) #type 10 dice si hay outlier

## 
##  Grubbs test for one outlier
## 
## data:  datos_prueba$W
## G = 5.86105, U = 0.96558, p-value = 1.7e-06
## alternative hypothesis: highest value 6 is an outlier

#Rechazamos H0, el maximo es outlier

grubbs.test(datos_prueba$W, type = 10, opposite = TRUE) #probamos con el minimo

## 
##  Grubbs test for one outlier
## 
## data:  datos_prueba$W
## G = 2.83693, U = 0.99194, p-value = 1
## alternative hypothesis: lowest value -2.89161896201262 is an outlier

#No rechazamos H0, el maximo no es outlier

grubbs.test(datos_prueba$W, type = 11) #type 11 dira si el menor y mayor son atipicos

## 
##  Grubbs test for two opposite outliers
## 
## data:  datos_prueba$W
## G = 8.69798, U = 0.95755, p-value = 0.0003314
## alternative hypothesis: -2.89161896201262 and 6 are outliers

#Rechazamos H0 el maximo y/o minimo son outliers

#Dixon
datos_prueba_dixon <- c(datos_prueba$W[1:23], max(datos_prueba$W))

dixon.test(datos_prueba_dixon)

## 
##  Dixon test for outliers
## 
## data:  datos_prueba_dixon
## Q = 0.57893, p-value < 2.2e-16
## alternative hypothesis: highest value 6 is an outlier

dixon.test(datos_prueba_dixon, opposite = TRUE)

## 
##  Dixon test for outliers
## 
## data:  datos_prueba_dixon
## Q = 0.18724, p-value = 0.9055
## alternative hypothesis: lowest value -2.22190001862843 is an outlier

#Rosner
rosnerTest(datos_prueba$Z, k = 10, alpha = 0.01)

## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 8.367234 8.576453 8.019786 7.447655 7.191904 5.923550 5.907384 3.385861 
##      R.9     R.10 
## 3.256228 3.216394 
## 
## $sample.size
## [1] 1000
## 
## $parameters
##  k 
## 10 
## 
## $alpha
## [1] 0.01
## 
## $crit.value
##  lambda.1  lambda.2  lambda.3  lambda.4  lambda.5  lambda.6  lambda.7  lambda.8 
##  4.396763  4.396529  4.396295  4.396061  4.395826  4.395592  4.395357  4.395122 
##  lambda.9 lambda.10 
##  4.394886  4.394650 
## 
## $n.outliers
## [1] 7
## 
## $alternative
## [1] "Up to 10 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##    [1]  3.455800e-01  1.266283e+00  1.554506e+00 -1.770120e+00  3.187243e-02
##    [6]  1.176828e+00 -8.714269e-01  1.118470e-01 -1.417057e+00  1.181479e+00
##   [11]  1.992825e+00 -2.345337e-01 -3.038342e-02  2.733995e+00  5.509904e-01
##   [16]  2.950506e-01 -7.256670e-01  5.800997e-01  7.955755e-01  1.002730e+00
##   [21] -8.397401e-01 -6.259758e-01  1.650710e+00  5.737262e-01  6.869517e-03
##   [26] -1.171198e+00 -4.674365e-01  7.384031e-01 -3.765558e-01 -1.474706e+00
##   [31] -5.913877e-01 -4.724967e-02 -4.798542e-01 -4.059122e-01 -1.469927e+00
##   [36] -6.267980e-01  1.062349e+00  2.095626e+00 -2.388879e-01 -7.151908e-02
##   [41]  6.717264e-01  4.679259e-01  5.118969e-01 -1.034128e-01  1.855991e+00
##   [46] -3.702029e-01 -5.599825e-01  4.623494e-02 -2.839046e-01  8.360382e-01
##   [51] -1.287253e+00  1.575307e-01 -7.249193e-01  2.653419e-01  1.661515e-01
##   [56] -5.824047e-01 -6.934094e-01 -1.535665e+00 -3.611391e-01  1.861446e+00
##   [61]  5.872332e-01  8.571039e-01 -1.574544e-02  7.437985e-01 -1.527491e+00
##   [66] -6.330369e-01  1.397439e+00  4.610322e-01 -8.170619e-01  8.229161e-01
##   [71] -9.531537e-02  2.031341e+00 -2.564079e-01 -1.579766e+00 -1.649734e-01
##   [76]  3.524168e-01 -7.574601e-01  8.737398e-04  5.811810e-04  2.081961e-01
##   [81] -3.450451e-01  5.302242e-01 -2.930053e-02 -7.288496e-01 -3.662936e-01
##   [86]  1.837857e-01 -1.344948e-01 -6.172187e-01  1.090672e-01 -1.537313e+00
##   [91]  7.248760e-01 -1.402322e+00 -5.512596e-01 -2.583508e+00  3.706332e-01
##   [96] -1.202274e+00  6.803633e-01 -4.441293e-01 -1.378279e+00 -1.311368e+00
##  [101]  6.407357e-01 -4.658736e-01  9.321953e-02  7.408514e-01  1.429357e+00
##  [106] -7.116064e-01 -2.921836e-01  9.752730e-01  1.555107e-01  2.835362e-01
##  [111] -3.325681e-01  1.982753e+00 -1.245377e-01 -1.206080e+00 -2.291880e-01
##  [116] -1.222117e+00  6.371631e-01  5.129341e-01 -5.691496e-01 -2.291711e+00
##  [121]  2.387171e-01 -2.336353e+00  4.817715e-01  4.738310e-01 -1.138586e+00
##  [126]  7.982135e-01  1.317077e+00  6.851203e-01  5.008756e-01 -1.350373e+00
##  [131]  1.055785e+00 -9.884559e-01 -2.197837e+00 -1.839086e+00  2.151893e-01
##  [136]  5.358055e-01 -1.314864e+00  2.765230e-01 -1.206791e+00  3.286068e+00
##  [141]  5.519513e-01 -1.892830e-01 -6.567023e-01 -2.834564e+00  7.231114e-02
##  [146]  3.854661e-01  2.308686e+00  6.406383e-01  4.234189e-01 -2.175787e+00
##  [151] -5.752590e-01 -1.612873e+00 -3.365808e-01  8.992614e-01  5.962451e-02
##  [156] -2.190877e+00  4.038998e-01  1.199579e+00 -4.837393e-01 -1.020747e+00
##  [161]  9.010070e-02 -1.079036e+00 -2.041726e-02  6.660272e-01  1.105933e-02
##  [166] -7.010801e-01  3.137635e+00 -1.305390e-01  8.816701e-01 -6.775299e-01
##  [171] -5.415201e-01  1.150277e-01  7.996112e-01  1.329789e-02  8.084004e-01
##  [176] -1.845699e-01 -3.152958e-01  1.266957e+00  5.922681e-01  1.202817e+00
##  [181]  2.582143e-01  5.494149e-02 -8.977975e-01 -8.243972e-01 -6.674073e-01
##  [186]  2.192055e-01 -4.221382e-01 -1.122868e+00  3.237151e-01 -2.394080e-02
##  [191]  7.687917e-01  1.283678e+00 -6.469983e-01 -4.634727e-02  5.229828e-01
##  [196]  6.274620e-01 -2.957732e-01 -9.008830e-01 -1.439175e-01 -8.957215e-01
##  [201]  2.248675e-01  3.097240e-01 -1.495337e+00  1.768583e-01 -2.574353e-01
##  [206]  2.022768e+00  7.149122e-01  3.370885e-01 -1.653952e+00  4.360470e-01
##  [211]  1.649348e+00  1.093541e+00  9.324829e-02  4.348040e-01  1.056065e+00
##  [216] -3.697729e-01  3.504336e-02 -9.474238e-01  2.793112e-01 -4.656577e-01
##  [221] -1.112773e+00 -8.568855e-01 -7.920995e-01  8.526145e-01  1.072185e+00
##  [226]  1.687825e+00  3.016684e-01 -6.766442e-01  6.380037e-01 -4.531747e-01
##  [231]  1.385558e+00  1.223366e+00 -1.697589e-01 -1.828745e+00  3.872882e-01
##  [236] -1.211891e+00 -1.034880e-01 -1.340877e-01  2.128408e+00  9.332909e-01
##  [241]  1.477739e+00 -7.985701e-01 -1.129282e+00 -6.335298e-01 -6.590816e-01
##  [246] -1.337504e-01  3.172361e-01 -1.612683e+00 -6.639627e-01 -9.371515e-01
##  [251]  1.194065e-01 -3.302462e-01 -1.323761e+00  7.878827e-01 -8.940506e-01
##  [256] -1.900288e+00 -9.170256e-01 -1.045441e+00 -9.241627e-01 -1.123301e+00
##  [261]  1.614284e-02  1.530006e+00  3.702515e-01 -8.835859e-02  3.894361e-02
##  [266]  1.545111e+00 -5.330771e-01 -6.791903e-01 -1.407241e+00 -6.617312e-01
##  [271]  6.829054e-02 -1.480176e-01  2.028525e+00  1.173808e-01  3.175251e-01
##  [276] -1.021428e+00  5.752402e-01  7.511665e-01 -1.099230e+00  1.816772e+00
##  [281]  5.100157e-01  2.396341e+00  1.743039e-01  1.109933e-01 -3.012849e-02
##  [286] -3.582338e-01 -9.470715e-01  3.892393e-01  2.119685e+00  9.225421e-01
##  [291]  9.475052e-01 -2.397436e-01 -2.036743e+00  8.277755e-02 -3.091358e-01
##  [296] -1.863835e+00 -1.547343e+00  1.841151e-01 -4.105728e-01  1.299636e+00
##  [301] -1.336311e+00 -4.110978e-01  6.153931e-01  1.174588e+00  4.492896e-01
##  [306]  1.698570e-01 -5.865420e-01 -1.830219e-01  6.630968e-01  9.577335e-01
##  [311]  1.275957e+00  1.653644e+00 -1.188735e+00  1.060806e+00 -1.213420e+00
##  [316]  4.857116e-02 -9.159486e-01  3.358809e-01 -1.060109e+00  1.972504e+00
##  [321] -8.762634e-01  3.784736e-01 -1.061032e+00 -1.785758e+00  3.067519e-01
##  [326] -1.309726e+00 -7.767471e-01  1.892410e-01  1.413067e+00  4.017924e-01
##  [331]  3.918412e-02  2.953814e-01 -5.165868e-01 -9.942932e-01 -6.935278e-01
##  [336]  6.879634e-01 -4.452412e-01  5.963613e-02  2.360411e-01 -1.101096e+00
##  [341] -5.960123e-01  6.327483e-01 -1.216377e+00 -2.395076e+00  1.295469e+00
##  [346]  1.387848e-01  1.493919e+00 -1.661893e-01 -2.522970e-01  2.031854e-01
##  [351] -4.280653e-01 -1.075752e+00  2.963791e-01 -9.429929e-01 -2.307329e-01
##  [356] -1.207815e+00  2.045574e-02  9.000102e-01 -3.326534e-01  8.925190e-01
##  [361] -6.373363e-01 -9.288169e-01 -4.143211e-01  1.294089e+00 -7.068576e-01
##  [366]  6.916990e-01  3.482572e-01 -1.049099e+00  5.740705e-01 -1.473523e+00
##  [371]  3.450773e-01 -7.751618e-02 -4.403495e-01  1.788593e+00  9.325019e-01
##  [376] -2.652675e-02  1.836423e+00  4.281765e-01  4.827890e-01  9.635802e-01
##  [381] -1.277553e+00 -1.767837e+00  9.084069e-01 -4.876830e-01 -5.210279e-01
##  [386]  1.119731e+00  3.918662e-01 -1.272538e+00 -9.245390e-01  4.925820e-01
##  [391]  4.130792e-01  3.237040e-01 -1.231492e-01 -3.092543e-01  2.030418e-02
##  [396]  2.369647e+00  5.918195e-01  1.171872e+00  3.001482e-01 -7.354368e-01
##  [401]  3.714252e-01  6.306462e-01  1.069918e+00 -1.871367e-01  1.075009e+00
##  [406] -2.216429e-01  6.379073e-01 -7.791417e-01 -6.972318e-01 -2.284589e-01
##  [411] -8.358460e-01  1.050842e+00 -1.749096e+00 -1.460876e+00  9.703787e-01
##  [416]  1.133196e+00  3.894688e-01  7.335964e-01 -1.362060e+00 -2.133380e+00
##  [421] -4.045231e-02 -5.102190e-02 -2.009850e+00  9.285540e-01 -3.965339e-01
##  [426] -1.886932e-01  3.403838e-01 -7.084488e-01  2.869804e-02 -1.380203e+00
##  [431] -8.056682e-01 -5.248658e-01  4.182708e-01  6.010931e-01 -8.997150e-03
##  [436] -4.194628e-02 -2.190535e+00 -8.448332e-01 -1.744446e-01  6.764618e-01
##  [441]  6.429242e-01  8.219183e-01  2.133424e-01 -1.204595e+00  9.021350e-01
##  [446] -5.981122e-01  3.449381e-01 -6.076412e-01 -3.942459e-01  1.454743e-01
##  [451] -9.412721e-01  3.463689e-01 -4.844844e-01 -1.101056e-02 -3.299023e-01
##  [456] -6.470943e-01  4.229783e-01 -6.554221e-01  1.494995e+00 -2.728105e+00
##  [461]  6.330653e-03 -1.338550e+00 -6.976849e-01 -1.321959e-01  6.119732e-01
##  [466]  1.671660e+00  4.431081e-01  1.328516e+00  1.796421e+00  4.005672e-02
##  [471] -1.297567e+00  1.720298e-01 -6.279012e-01  2.053544e-01 -2.757870e-01
##  [476]  4.937983e-01 -6.842283e-03  7.283671e-01  5.455180e-01 -2.102092e+00
##  [481] -8.508022e-01 -8.938886e-01  1.220790e+00  1.467932e-01 -4.903006e-01
##  [486]  7.451385e-01  1.164427e+00 -3.328572e-01 -1.048796e-03  2.324123e-01
##  [491]  2.448543e+00  5.071383e-01  6.278285e-01  4.364941e-01  4.367971e-01
##  [496] -7.733525e-01 -4.712824e-01 -9.973433e-01  1.856334e+00 -1.959507e+00
##  [501] -1.967565e+00 -6.824346e-01 -6.493233e-02 -6.499535e-02  1.383887e+00
##  [506] -6.105673e-01  1.285755e+00 -9.715797e-01  1.352246e-01  6.975456e-01
##  [511] -1.793490e+00 -1.321117e-01 -6.357379e-01 -1.777088e+00  1.825485e-02
##  [516] -5.930668e-01  2.391676e-02 -2.126900e+00 -3.005899e+00  5.927737e-01
##  [521]  9.675236e-02 -9.091856e-02 -1.428364e-01  2.187257e-01 -9.222978e-01
##  [526] -4.281314e-01  1.331583e+00  1.387693e+00 -1.690462e+00 -5.253478e-01
##  [531]  7.741131e-01 -3.322174e-01  7.231831e-01  1.488301e-01 -1.924902e+00
##  [536]  1.559708e-01  2.127099e+00 -1.839418e+00  6.512195e-02 -2.316673e-01
##  [541]  1.549987e-01 -3.933221e-01  3.397903e-01  4.153700e-01 -8.918293e-01
##  [546]  8.441103e-01  1.358814e+00 -1.563920e+00  1.965194e-01 -6.623831e-01
##  [551] -2.634806e+00  2.088810e-01  1.613975e+00 -1.501161e+00 -1.546683e+00
##  [556]  1.394478e+00  3.310006e-02  3.978664e-02  6.235650e-01  1.397165e+00
##  [561]  2.246111e-01  1.209143e+00  6.216189e-02  9.536513e-01  1.098242e+00
##  [566] -1.867337e+00 -1.093118e+00 -5.881632e-01 -1.669563e+00 -4.221590e-01
##  [571]  1.000239e+00 -1.675040e+00 -6.197989e-01  6.009983e-01 -1.941745e+00
##  [576] -4.864742e-01 -1.654852e+00  3.260689e-01 -2.018057e+00 -5.618141e-01
##  [581]  7.192486e-01 -6.427769e-01 -1.337600e+00  2.221531e-01  1.601590e+00
##  [586]  4.318319e-01  2.526529e-01 -3.405555e-01 -7.787286e-01  2.971440e-01
##  [591]  9.184401e-01 -8.124326e-02 -2.639220e-02  8.044779e-01 -9.399841e-01
##  [596] -2.097389e+00 -9.305081e-01  1.419864e-01 -5.818625e-01 -8.841333e-01
##  [601]  1.355435e+00 -5.505211e-02 -1.215141e+00  8.626272e-01  4.858170e-01
##  [606] -2.839180e-01  2.071145e-01 -6.455851e-01  6.630249e-01 -3.550153e-01
##  [611]  1.893842e-02  1.217350e+00 -9.780026e-01  1.939094e+00 -7.189517e-01
##  [616] -5.373968e-02  7.392162e-01 -1.733626e-01  2.601406e-02 -7.120083e-01
##  [621] -7.639327e-01 -3.267842e-01 -9.422777e-02 -1.255360e+00 -6.559904e-01
##  [626]  2.791439e-01  8.171340e-01  5.791277e-01  3.510311e-01 -1.491407e-01
##  [631]  1.365120e+00  1.398977e-02  1.309718e+00 -8.338261e-02  1.579668e-01
##  [636]  2.085417e-01  3.024028e-01  8.702331e-01 -3.541082e-01 -4.047706e-01
##  [641] -1.344618e+00 -1.262980e+00  1.867794e+00  3.374351e-01 -1.600216e+00
##  [646] -2.091152e-01 -3.424309e-01  1.281903e+00 -9.789756e-01 -9.620664e-01
##  [651] -2.297359e+00 -9.067311e-01 -1.764399e+00 -1.346051e+00  7.585429e-01
##  [656]  3.606275e-01 -3.872997e-01 -4.402435e-01  1.009058e-01 -4.644074e-01
##  [661]  8.548708e-01  2.960614e-01 -3.413563e-01  1.597808e-01  6.419324e-01
##  [666]  1.434272e-01  1.646396e+00 -1.349010e-01  1.443222e+00 -8.468319e-01
##  [671]  1.458524e+00  4.188131e-01  4.709365e-01 -3.279470e-01 -2.529762e-01
##  [676]  9.208755e-01  5.374249e-02 -8.400624e-01 -8.270101e-01 -8.924386e-01
##  [681]  2.144497e+00 -1.030984e-01 -1.112351e+00 -3.102044e-01  6.289800e-01
##  [686]  1.314286e+00 -2.367430e-01  5.437108e-01 -1.254036e+00  1.248800e-01
##  [691]  7.240603e-01  1.538292e+00 -1.048794e+00 -5.967742e-01 -8.265588e-02
##  [696]  5.505766e-01  6.980699e-02 -1.890128e+00 -1.913347e-01 -1.816804e+00
##  [701] -1.022407e+00  2.371683e+00 -6.849838e-01  2.079055e-01  2.830204e-01
##  [706]  1.546155e-01  1.451134e-01  8.131586e-01  4.809818e-01 -1.887754e+00
##  [711] -1.051132e+00 -1.385061e+00 -8.135246e-01 -2.317169e+00 -1.904642e+00
##  [716]  9.958018e-01  6.996124e-01  1.320830e+00  1.667857e+00 -3.713937e-01
##  [721]  1.622926e-01 -3.422973e-01 -1.222366e+00  1.273161e-01 -2.766362e-01
##  [726] -5.842587e-01 -9.827966e-01 -6.565223e-02  6.291359e-01  5.760232e-01
##  [731]  1.057040e+00  2.324090e-01  8.762689e-01 -6.130374e-01  4.080992e-02
##  [736]  9.470359e-01  4.089747e-01 -8.480561e-01  4.396661e-01  6.556652e-01
##  [741] -2.261269e-01 -1.402236e+00 -9.473668e-01  7.143116e-01  1.089909e+00
##  [746]  5.721682e-01 -1.146955e+00 -1.778124e-01 -4.293113e-01 -8.776319e-01
##  [751] -4.973774e-01  4.279697e-01 -3.125153e-01  7.071850e-01 -2.637680e-01
##  [756]  3.197695e-01 -1.422258e-01 -2.024245e+00 -6.871062e-01 -9.557765e-02
##  [761]  1.659793e+00 -1.748481e+00  5.655623e-01 -9.363930e-01 -2.809868e-01
##  [766] -1.877386e-01  7.022130e-01  6.189342e-01  7.869553e-01 -5.879198e-01
##  [771] -2.723514e+00 -1.269503e-02  2.081825e+00 -7.077014e-01  2.836899e-01
##  [776] -1.286041e+00  2.203286e-01 -1.305333e+00  6.606951e-01 -9.939146e-01
##  [781] -3.269100e-01 -2.835109e-01  8.592673e-01  5.406038e-01  2.016749e+00
##  [786] -1.590069e+00  1.507455e+00  7.973404e-02 -1.040812e+00  1.522838e-01
##  [791] -4.369705e-01  3.545981e-01  1.198017e+00 -4.686640e-01  4.713734e-01
##  [796] -4.116314e-01  1.454336e+00  2.175627e+00  6.079554e-01 -1.168654e+00
##  [801]  1.057811e+00  1.266839e-01  2.877983e-01 -9.458880e-01 -9.957209e-01
##  [806]  6.555010e-01 -4.459709e-02 -1.003693e+00 -1.385690e+00  1.337810e-01
##  [811] -3.692614e-01  1.906247e+00 -4.988397e-01  2.067150e+00  6.440565e-01
##  [816] -6.504902e-01 -8.541471e-02 -1.143555e+00 -8.420595e-01 -1.167957e+00
##  [821]  1.272502e+00 -4.823582e-01 -4.270786e-01  2.459491e-01  8.212876e-01
##  [826] -1.332622e+00  1.519973e+00 -6.207134e-01  1.673938e-01 -9.321678e-01
##  [831]  6.392905e-01 -8.298971e-01  1.655125e+00 -2.254584e-01  2.597423e+00
##  [836] -1.249643e+00 -6.685463e-01 -7.392952e-01  6.530875e-01 -1.648756e+00
##  [841]  6.288590e-01 -1.381345e-01  2.341478e-01 -1.142839e-01  1.189733e+00
##  [846]  8.557177e-01 -4.744168e-01 -1.146835e+00  2.643036e-01 -1.220340e+00
##  [851] -1.488070e+00  2.289226e+00  1.224135e+00  2.762064e-01 -1.101104e+00
##  [856]  3.260150e-01 -6.188469e-01 -5.744232e-01  5.919254e-02 -4.039209e-02
##  [861] -7.229437e-01  6.109809e-01  9.128689e-01 -3.243701e-01 -3.726119e-01
##  [866]  2.147467e-02  6.632164e-01 -3.184063e-01 -1.167085e+00 -1.123876e-01
##  [871]  6.731182e-01  1.608554e+00 -1.172903e+00 -7.455932e-03 -3.056263e-01
##  [876] -2.110416e+00 -6.253476e-01 -5.897900e-01  5.972193e-02 -2.155201e+00
##  [881] -6.993075e-01  5.511627e-01  4.889365e-02 -8.975387e-01  1.774632e-01
##  [886]  1.421493e-01  1.026265e+00 -2.631159e-01 -3.208471e+00  2.893501e-01
##  [891] -6.873571e-01  6.767711e-01 -2.871863e-02  6.551471e-01 -1.046066e+00
##  [896]  7.941311e-01  4.550781e-01 -4.788560e-01 -1.094614e+00 -7.570619e-01
##  [901]  1.468447e+00 -1.522762e+00  6.731775e-01  1.817797e+00  2.020632e-01
##  [906] -1.199253e+00  9.462436e-01 -1.200680e-01  9.850433e-01  1.413268e-01
##  [911] -2.339021e-01 -1.321681e+00  8.500955e-01  1.227687e-01  7.664268e-01
##  [916] -5.782485e-01  4.045617e-01  3.065375e-01 -3.503417e-01  1.542357e+00
##  [921]  5.806893e-01 -5.715692e-01 -7.372512e-01  1.094498e-01  6.887060e-02
##  [926] -1.045290e+00  1.398563e+00  9.043678e-01  2.432016e-01  8.445182e-01
##  [931]  1.192510e+00 -1.602163e+00 -7.341609e-01  3.476621e-01  1.900018e+00
##  [936] -6.608515e-01 -7.911563e-01 -1.758122e+00 -1.933416e+00 -1.535481e+00
##  [941]  2.379056e-01 -5.675204e-01  1.657937e+00 -8.294604e-01  5.617480e-01
##  [946]  1.656227e+00 -1.003164e+00 -1.686554e+00 -2.137754e-01  4.510996e-01
##  [951] -9.473775e-01  1.039101e-01  4.265191e-01 -3.771758e-01  5.660057e-01
##  [956] -5.571301e-01 -4.855463e-01 -5.233411e-01  2.941783e-01  1.040821e-01
##  [961] -1.791602e-01  6.390583e-02  1.883583e-01  1.220886e+00  1.251114e+00
##  [966]  3.394930e-01 -1.538456e+00  3.546281e-01 -8.399735e-01 -8.784880e-01
##  [971] -6.620407e-01  2.123615e+00  4.071828e-01 -1.358478e+00 -3.341674e-01
##  [976] -1.821440e+00  2.296234e-01  1.646561e-01 -5.899871e-01 -6.400729e-01
##  [981]  1.394082e+00  1.552965e-01  1.030991e+00 -9.411904e-01  5.467628e-01
##  [986] -6.994582e-01  6.949818e-02 -1.626145e-01  9.586295e-02  1.308300e-01
##  [991] -9.357874e-01  1.438972e+00 -1.043152e+00 -9.000000e+00 -6.000000e+00
##  [996] -1.000000e+01  1.000000e+01  8.000000e+00  6.000000e+00  7.500000e+00
## 
## $data.name
## [1] "datos_prueba$Z"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##    i      Mean.i      SD.i      Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1  0 -0.05078390 1.2012075  10.000000     997 8.367234   4.396763    TRUE
## 2  1 -0.06084475 1.1588888 -10.000000     996 8.576453   4.396529    TRUE
## 3  2 -0.05088567 1.1158795  -9.000000     994 8.019786   4.396295    TRUE
## 4  3 -0.04190963 1.0797908   8.000000     998 7.447655   4.396061    TRUE
## 5  4 -0.04998384 1.0497893   7.500000    1000 7.191904   4.395826    TRUE
## 6  5 -0.05757176 1.0226252   6.000000     999 5.923550   4.395592    TRUE
## 7  6 -0.06366590 1.0049007  -6.000000     995 5.907384   4.395357    TRUE
## 8  7 -0.05768772 0.9875643   3.286068     140 3.385861   4.395122   FALSE
## 9  8 -0.06105844 0.9823308   3.137635     167 3.256228   4.394886   FALSE
## 10 9 -0.06428618 0.9775495  -3.208471     889 3.216394   4.394650   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"

#Probamos la existencia de 10 outliers, pero como pusimos 7 manualmente, dice que 
#3 son falsos

Ahora, una aplicación de la detección de anomalías es la ley de Benford, que mencionamos previamente. Esta ley se utiliza para la detección de fraudes, un buen modelo de detección de fraudes debe cumplir con:

Buen accuracy
Fácil de interpretar
Cumplimiento normativo
Costo razonable
Complementar los enfoques basados en expertos

Volviendo a la ley de Benford, esta es una de las herramientas más sencillas y más utilizadas para la detección de fraudes, y es conocida como la ley de los primeros dígitos. Esta ley establece un patrón que el primer dígito de una serie de valores debe seguir, este primer dígito debe ser el dígito más a la izquierda. Mientras que uno piensa que la distribución es uniforme, la ley de Benford espera que la distribución sea esta:

Para que la ley se cumpla se deben seguir las siguientes condiciones:

Se examinara una variable numérica positiva que no puede tomar el valor de 0
Las observaciones son generadas aleatoriamente o provienen de un proceso que tiene cierto grado de aleateoriedad.
Los valores no se encuentran restringidos por un máximo o un mínimo
La muestra es grande
El orden de magnitud debe ser relativamente grande

Con esto dicho, la ley de Benford se aplica analizando el orden de magnitud y orden de magnitud robustos si ambos son mayores a 3, o uno es mayor a 3 y el otro esta cerca a 3 se puede aplicar:

#Evaluación preliminar
OOM <- log10(max(datos6$Amount)/min(datos6$Amount)) 
OOM #orden de magnitud

## [1] 7.427543

ROM <- log10(quantile(datos6$Amount, 0.99)/quantile(datos6$Amount, 0.01))
ROM #orden de magnitud robusto

##      99% 
## 3.702056

#ambos valores son superiores a 3, lo que significa que tiene un orden de magnitud suficiente
#para aplicar la ley de benford

## Analisis Grafico
ben <- benford(datos6$Amount, number.of.digits = 1)

plot(ben, except = c("second order", "summation", "mantissa", "chi squared", 
                     "abs diff", "ex summation", "Legend"), multiple = FALSE)

suspectsTable(ben)

##    digits absolute.diff
##     <int>         <num>
## 1:      5      4093.028
## 2:      9      3710.770
## 3:      1      3245.456
## 4:      2      2813.341
## 5:      3      2756.202
## 6:      4      2602.369
## 7:      7      1500.411
## 8:      6      1241.791
## 9:      8       135.139

Con esto realizamos preparamos todo para el análisis estadístico, para este debemos revisar si hay una diferencia significativa entre lo real y lo esperado, es decir, los datos que se tienen y la ley de Benford, si lo hay eso puede significar un posible fraude:

Para el análisis podemos realizar 3 cosas:

Análisis gráfico de Intervalos de Confianza para los dígitos: si los datos quedan dentro del intervalo que el algoritmo propone, se pueden concluir que esos dígitos siguen la ley de Benford

signifd.analysis(datos6$Amount, digits = 1, ci_lines = c(0.01))

## $summary
##               1            2            3            4         5            6
## freq  0.3186289 1.608355e-01 1.099928e-01 8.279830e-02 0.1013763 6.021300e-02
## pvals 0.0000000 2.671834e-66 6.225897e-84 2.962559e-93 0.0000000 5.843537e-31
##                  7          8          9
## freq  4.985576e-02 0.05041971 0.06587966
## pvals 1.587332e-50 0.15317256 0.00000000
## 
## $CIs
##               1         2         3          4          5          6          7
## 0.005 0.2982786 0.1738065 0.1229554 0.09513553 0.07756160 0.06544766 0.05658999
## 0.5   0.3010300 0.1760913 0.1249387 0.09691001 0.07918125 0.06694679 0.05799195
## 0.995 0.3037814 0.1783760 0.1269220 0.09868450 0.08080089 0.06844592 0.05939390
##                8          9
## 0.005 0.04983106 0.04450411
## 0.5   0.05115252 0.04575749
## 0.995 0.05247398 0.04701087

Chi Cuadrado: se desea ver si los datos reales siguen lo esperado, para esto hay dos hipótesis, la H0 dice que los datos se comportan como la ley de Benford y una HA que dice que los datos no se comportan según lo esperado

chisq(ben)

## 
##  Pearson's Chi-squared test
## 
## data:  datos6$Amount
## X-squared = 4258.6, df = 8, p-value < 2.2e-16

En este caso los datos no se comportan según lo esperado (rechazamos H0)

Analisis MAD: nos dirá que tan fuerte es la relación entre los datos y la ley de Benford, siguiendo estas interpretaciones.
- Relación estrecha: entre 0 y 0.006
- Relación aceptable: entre 0.006 y 0.012
- La relación es marginal: entre 0.012 y 0.015
- No hay relación: mayor a 0.015

MAD(ben)

## [1] 0.0133147

En este caso, la relación es marginal. Por ultimo se pueden tomar los números más fraudulentos que en este caso son 5 y 9 y ver las observaciones:

ben_5_9 <- getDigits(ben, datos6, digits = c(5, 9))
head(ben_5_9)

##    registro VendorNum    Date Amount
##       <int>     <int>  <char>  <num>
## 1:        4      2001 2/01/10  59.00
## 2:        5      2001 2/01/10  59.56
## 3:        6      2001 2/01/10  50.38
## 4:       12      2001 2/01/10  55.22
## 5:       15      2001 2/01/10  55.29
## 6:       19      2001 2/01/10  94.29

Y se pueden ver los duplicados:

duplicatesTable(ben)

##         number duplicates
##          <num>      <int>
##     1:   50.00       6022
##     2: 1153.35       2264
##     3: 1083.45       1185
##     4:  150.00       1056
##     5:  988.35       1018
##    ---                   
## 65425:   38.22          1
## 65426:    5.36          1
## 65427:   28.49          1
## 65428:  424.28          1
## 65429:  108.97          1

duplicados <- getDuplicates(ben, datos6, how.many = 2)
duplicados

##       registro VendorNum     Date Amount
##          <int>     <int>   <char>  <num>
##    1:     8930      2098 31/08/10     50
##    2:    10963      2188  9/12/10     50
##    3:    11310      2208 17/09/10     50
##    4:    11312      2209 20/01/10     50
##    5:    11846      2245 23/04/10     50
##   ---                                   
## 8282:   176565     21798  4/02/10     50
## 8283:   176585     21814 22/02/10     50
## 8284:   176781     21960 30/01/10     50
## 8285:   180578     34061 14/02/10     50
## 8286:   180600     34083 14/02/10     50

Repaso Parcial Final